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PREFACE 


Many important problems of complexity are related in one way or 
another with the presence of phase transition phenomena. Most 
complex systems are known to potentially display a number of 
different patterns of qualitative behavior or phases. Such phases 
correspond to different forms of internal organization and two 
given phases are usually separated by a sharp boundary, and 
crossing such a frontier implies a change in system-level behavior. 
Many of these transitions reveal important features of the system 
undergoing them. A good model of a complex system should be 
able to predict the presence of such phases and their implications. 
Some of these changes can be catastrophic, and understanding 
them is crucial for the future of biodiversity or even our society. 
The existence of different phases that the same system can 
display has been well known in physics for centuries. That 
water freezes and becomes solid or boils to become a gas are 
part of our daily experience with matter. And yet, a close 
examination of these events immediately triggers deep questions: 
Why are these changes sharp? Why not simply a continuous 
transformation? What drives these changes at the microscopic 
level? How much of these details are needed to understand the 
phenomenon at the large scale? The presence of phase transitions 
in ecology, economy, epidemics, or cancer, just to cite a few 


xii Preface 


examples discussed in this book, has been a matter of analysis 
and theoretical work over the last decades. Examples of them 
are scattered over the technical literature and I have tried to put 
together a relevant group of examples. 

This book is not a general treatise of phase transitions, which 
would require a rather thick volume, given the multiple aspects 
required to cover the field. Instead, I have chosen a particular 
approximation, —the so-called mean field theory—and in partic- 
ular those cases where such theory allows simplification into the 
most basic mathematical model. Mean field theory ignores some 
crucial aspects of reality, such as the local nature of interactions 
among elements or the role played by fluctuations. In removing 
these components we can formulate the simplest possible model. 

Many people have contributed to this book. I would like 
to thank colleagues that have crossed my path during this long 
journey, particularly those at the Santa Fe Institute, where much 
of this book was written. I am also indebted to the members of 
the Complex Systems Lab for stimulating discussions. The book 
contains a list of references to help the reader further explore the 
problems presented here and that are beyond the limited scope 
of the mean field picture. The list is not exhaustive but contains 
books and papers that helped me understand the problems and 
look beyond them. I hope this book becomes part of that list. 
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This page intentionally left blank 


1 


PHASE CHANGES 


1.1 Complexity 


By tracking the history of life and society, we find much 
evidence of deep changes in life forms, ecosystems, and civiliza- 
tions. Human history is marked by crucial events such as the 
discovery of the New World in 1492, which marked a large- 
scale transformation of earth’s ecology, economics, and culture 
(Fernandez-Armesto 2009). Within the context of biological 
change, externally driven events such as asteroid impacts have also 
triggered ecosystem-scale changes that deeply modified the course 
of evolution. One might easily conclude from these examples that 
deep qualitative changes are always associated with unexpected, 
rare events. However, such intuition might be wrong. Take 
for example what happened around six thousand years ago in 
the north of Africa, where the largest desert on our planet 
is now located: the Sahara. At that time, this area was wet, 
covered in vegetation and rivers, and large mammals inhabited 
the region. Human settlements emerged and developed. There 
are multiple remains of that so-called Green Sahara, including 
fossil bones and river beds. The process of desertification was 
initially slow, when retreating rains changed the local climate. 
However, although such changes were gradual, at some point 
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the ecosystem collapsed quickly. The green Sahara became a 
desert. 

Transitions between alternate states have been described in 
the context of ecology (Scheffer 2009) and also in other types 
of systems, including social ones. Complex systems all display 
these types of phenomena (at least as potential scenarios). But 
transitions can also affect, sometimes dramatically, molecular 
patterns of gene activity within cells, behavioral patterns of 
collective exploration in ants, and the success or failure of cancer 
or epidemics to propagate (Solé et al. 1996; Solé and Goodwin 
2001). When a given parameter is tuned and crosses a threshold, 
we observe change in a system’s organization or dynamics. We 
will refer to these different patterns of organization as phases. The 
study of complexity is, to a large extent, a search for the principles 
pervading self-organized, emergent phenomena and defining its 
potential phases (Anderson 1972; Haken 1977; Nicolis and 
Prigogine 1977, 1989; Casti 1992a, b; Kauffman 1993; Cowan 
et al. 1994; Gell-Mann 1994; Coveney and Highfield 1995; 
Holland 1998; Solé and Goodwin 2001; Vicsek 2001; Mikhailov 
and Calenbuhr 2002; Morowitz 2002; Sornette 2004; Mitchell 
2009). Such transition phenomena are collective by nature and 
result from interactions taking place among many interacting 
units. These can be proteins, neurons, species, or computers (to 
name just a few). 

In physics, phase changes are often tied to changes between 
order and disorder as temperature is tuned (Stanley 1975; Binney 
et al. 1992; Chaikin and Lubensky 2001). Such phase transitions 
typically imply the existence of a change in the internal symmetry 
of the components and are defined among the three basic types 
of phases shown in figure 1.1. An example of such transition 
takes place between a fluid state, either liquid or gas, and a 
crystalline solid. The first phase deals with randomly arranged 
atoms, and all points inside the liquid or the gas display the same 
properties. In a regular (crystalline) solid, atoms are placed in the 
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G 


Figure 1.1. Three standard phases of matter: solid (S), liquid (L) 
and gas (G). Increasing a given external parameter (the control 
parameter) such as temperature, we can increase the degree of 
disorder. Such disorder makes molecules fluctuate around their 
equilibrium positions, move around them, or just wander freely. 


nodes of a regular lattice. In the gas phase (at high temperature) 
kinetic energy dominates the movement of particles and the 
resulting state is homogeneous and isotropic. All points are 
equivalent, the density is uniform, and there are essentially no 
correlations among molecules. In the liquid phase, although still 
homogeneous, short-distance interactions between molecules 
leads to short-range correlations and a higher density. Density 
is actually the fundamental difference distinguishing these two 
phases. Finally, the ordered arrangement observed at the solid 
phase is clearly different in terms of pure geometry. Molecules 
are now distributed in a highly regular way. Crystals are much 
less homogeneous than a liquid and thus exhibit less symmetry. 
However, beyond the standard examples of thermodynamic 
transitions between these three phases, there is a whole universe. 
Matter, in particular, can be organized in multiple fashions, and 
this is specially true when dealing with so-called soft matter. 
But the existence of different qualitative forms of macroscopic 
organization can be observed in very different contexts. Many are 
far removed from the standard examples of physics and chemistry. 


And yet, some commonalities arise. 


4 Phase Changes 
1.2 Phase Diagrams 


Phase changes are well known in physics: boiling or freezing 
are just two examples of changes of phase where the nature of 
the basic components is not changed. The standard approach 
of thermodynamics explores these changes by defining (when 
possible) the so-called equation of state (Fermi 1953), which is a 
mathematical expression describing the existing relations among a 
set of state variables (i.e., variables defining the state of a system). 
For an ideal gas, when no interactions among molecules need be 
taken into account,' the equation reads 


pV=nRT (1.1) 


where p, V, and T indicate pressure, volume, and temperature, 
respectively. Here R is the so-called gas constant (the same for 
all gases) and n the amount of substance (in number of moles). 
This equation is valid for a pure substance, and as we can see, it 
establishes a well-defined mathematical relation between p, V, T, 
and m. Given this equation, only three independent variables 
are at work (since the fourth is directly determined through the 
state equation). From this expression, we can plot, for example, 
pressure as a function of V and 7: 


nRT 
p(V, nar (1.2) 


which describes a continuous surface, displayed in figure 1.2a. 
For a given amount of perfect gas n, each point on this surface 
defines the only possible states that can be observed. Using this 
picture, we can consider several special situations where a given 
variable, such as temperature, is fixed while the other two can be 


! In this approximation, the molecules are considered point objects. However, 
in spite of this extreme oversimplification, the model can be used in a number 
of applied contexts. 


Figure 1.2. The equation of state describes the relations between the 
three key variables (p, V, T) that allow definition of a surface. In (a) 
the corresponding surface p(V, T) for an ideal gas is shown. This 
is obtained from the equation of state pV = nRT, which gives a 
continuous surface associated to a single gas phase. In reality (b) 
things are much more complex and the surface is discontinuous. This 
second plot corresponds to a substance (such as water) that expands 
when it solidifies. At high temperatures it is steam. Here kinetic 
energy is larger than the potential energy. But once T is decreased, 
nonlinear changes occur, as is evident from the discontinuous shape 
of the surface. 


changed. Fixing a given temperature Tj, we can see, for example, 
that pressure decays with volume following an inverse law, that 
is, p(V; Ti) =2RT\/V, which allows defining an isothermal 
process as the one moving on this line. The curve is called 
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an isotherm. By using different temperatures, we can generate 
different isotherms (which are in fact cross sections of the previous 
surface). Similarly, we can fix the volume and define another 
set of curves, now given by p(7) = BT with B = nR/V. The 
important idea here is that all possible states are defined by the 
equation of state and that in this case all possible changes are 
continuous. 

The previous equations are valid when we consider a very 
diluted gas at high temperature. However, in the real world, tran- 
sitions between different macroscopic patterns of organization 
can emerge out of molecular interactions” (figure 1.2b). Different 
phases are associated with different types of internal order and 
for example when temperature is lowered, systems become more 
ordered. Such increasing order results from a dominance of 
cohesion forces over thermal motion. Different combinations of 
temperature and pressure are compatible with different phases. 
An example is shown in figure 1.3, where the phase diagram for 
water is displayed? involving liquid, solid, and gas phases. Two 
given phases are separated by a curve, and crossing one of these 
curves implies a sharp change in the properties of the system. 
We say that a first order transition occurs and, in this case, each 
phase involves a given type of organization though no coexistence 
between phases is allowed. The melting of a solid (such as ice) or 
the boiling of a liquid (such as water) is a daily example, where 
both solid and liquid are present at the melting point. 

There is a special point in the previous diagram that appears 
when we follow the liquid-gas boundary curve and defines its 
limit. This so-called critical point describes a situation where there 
is in fact no distinction between the two phases. Moreover, we 


? This diagram is in fact just a small subset of the whole spectrum of possible 
phases displayed by water, which exhibits many other exotic phases (Stanley et al. 
1991; Ball 1999). 

3 Let us notice that this is not a completely typical phase diagram for a liquid, 
since the fusion curve has positive slope. 
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Figure 1.3. Phase diagram for water. Three basic phases are defined, 
separated by well-defined critical curves. A special path is indicated 
around the critical point that allows moving from one phase to 
the other (here from gas to liquid) without crossing any critical 
curve. Figures (b) and (c) illustrate the two different types of phase 
transitions that can occur. 


can see that the lack of a boundary beyond the critical point 
makes possible continuous movement from one phase to the other, 
provided that we follow the appropriate path (such as the one 
indicated in figure 1.3 with a dashed curve). The presence of 
this point has a crucial relevance in understanding the nature and 
dynamics of many natural and social phenomena. 


1.3 Interactions Make a Difference 


Phase transitions have been shown to occur in many different 
contexts (Chaikin and Lubensky 1995; Stanley et al. 2000; 
Solé and Goodwin 2001). These include physical, chemical, 
biological, and even social systems. In figure 1.4 we illustrate 
several examples, including those from molecular, behavioral, and 
cellular biology and cognitive studies. These are systems spanning 
many orders of magnitude in their spatial embodiment, and the 


Figure 1.4. Phase transitions can be identified in many differ- 
ent contexts. So (a-b), proteins, for example, starting from a 
linear chain (a) fold into three-dimensional objects (b) able to 
perform given functions. Amphiphilic molecules (such as lipids) 
spontaneously form complex spatial structures in water. De- 
pending on the concentration of molecules, different patterns 
are obtained (c) including layers and closed vesicles (image 
from http://en.wikipedia.org/wiki/Micelle). Swarm behavior in fish 
schools (d) can be modeled by means of different types of 
agent-based approaches. By changing the right parameters, we 
can see transitions between different types of collective motion 

(continued) 
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nature of the transitions is very different in terms of its functional 
and evolutionary relevance, but all of them have been described 
by means of simple models. 

Proteins exhibit two basic classes of spatial organization. Either 
they are folded in three-dimensional space or remain unfolded 
(as a more or less linear chain). The change from one state to 
the other (a—b) takes place under critical conditions (temperature 
or even the presence of other molecules). Our second example 
deals with a special and important group of molecules, so-called 
amphiphiles, which possess both hydrophilic and hydropho- 
bic properties. Typically, there is a hydrophilic head part (see 
figure 1.4c) that gets in touch with water molecules, whereas the 
opposing side, the so-called hydrophobic tail, avoids interacting 
with water. As a consequence, sets of amphiphilic molecules will 
tend to form self-organized, spatial structures that minimize the 
energy of interaction (Evans and Wennerstrém 1999; Mouritsen 
2005). By depending on the relative concentrations of water 
and amphiphiles, different stable configurations will be observed, 
defining well-defined phases (such as layers or vesicles). 

Our third and fourth examples involve behavioral patterns of 
interactions among individuals within groups of animals, as these 
occur with a fish school (figure 1.4d) or, at a lower scale of organi- 
zation, among cells in a cell culture (figure 1.4e). In the first case, 
interactions involve repulsion, attraction, the tendency to main- 
tain the same movement as neighboring individuals or to remain 
isolated. This is easily modeled using computer simulations and 
mathematical models. By tuning the appropriate parameters, we 


Figure 1.4 (continued). Bacterial colonies can display markedly 
different spatial patterns of growth (e) when the concentration of 
a given nutrient is tuned (modified from Cohen et al., 1996). And 
finally, (f) phase transitions occur in our brains as we shift from one 
image (a vase) to the alternate one (two faces) while looking at an 
ambiguous picture. (Pictures courtesy of Jacques Gautrais) 


10 Phase Changes 


can observe different types of collective motion and colony- 
level patterns (Vicsek et al. 1995; Toner and Tu 1998; Czirok 
et al. 1999; Camazine et al. 2001; Mikhailov and Calenbuhr 
2002; Theraulaz et al. 2003; Schweitzer 2003; Sumpter 2006; 
Garnier et al. 2007; Gautrais et al. 2008). Growing bacterial 
cell populations in a Petri dish with variable concentrations of 
a critical nutrient also demonstrate that different forms of colony 
organization emerge once given thresholds are crossed (Ben Jacob 
2003; Ben Jacob et al. 2004; Kawasaki et al. 1997; Matushita 
et al. 1998; Eiha et al. 2002; Wakano et al. 2003). The final case 
(figure 1.4f) involves cognitive responses to ambiguous figures 
(Attneave 1971) and a dynamical example of transitions among 
alternate states (Ditzinger and Haken 1989). Our brain detects 
one of the figures (the two faces) or the other (the vase), and both 
compete for the brain’s attention, which alternates between the 
two “phases” (Kleinschmidt et al. 1998). 

In this book we present several examples of phase transitions 
in very different systems, from genes and ecosystems to insect 
colonies or societies. Many other examples can be mentioned, 
including chemical instabilities (Nicolis et al. 1976; Feinn and 
Ortoleva 1977; Turner 1977; Nitzan 1978), growing surfaces 
(Barabdsi and Stanley 1995), brain dynamics (Fuchs et al. 1992; 
Jirsa et al. 1994; Kelso 1995; Haken 1996, 2002, 2006; Steyn- 
Ross and Steyn-Ross 2010), heart rate change (Kiyono et al. 
2005), immunology (Perelson 1989; de Boer 1989; Tomé and 
Drugowich 1996; Perelson and Weisbuch 1997; Segel 1998), 
galaxy formation (Schulman and Seiden 1981, 1986) cosmologi- 
cal evolution (Guth 1999; Linde 1994), computation (Landauer 
1961; Huberman and Hogg 1987; Langton 1990; Monasson 
et al. 1999; Moore and Mertens 2009), language acquisition 
(Corominas-Murtra et al. 2010), evolution of genetic codes 
(Tlusty 2007), politics and opinion formation (Lewenstein et al. 
1992; Kacperski and Holyst 1996; Weidlich 2000; Schweitzer 
2002, 2003; Buchanan 2007), game theory (Szabo and Hauert 
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2003; Helbing and Yu 2009; Helbing and Lozano 2010), and 
economic behavior (Krugman 1996; Arthur 1994, 1997; Ball, 
2008; Haldane and May 2011) to cite just a few. 


1.4 The Ising Model: From Micro to Macro 


In the previous sections we used the term critical point to describe 
the presence of a very narrow transition domain separating 
two well-defined phases, which are characterized by distinct 
macroscopic properties that are ultimately linked to changes in 
the nature of microscopic interactions among the basic units. A 
critical phase transition is characterized by some order parameter 
(mu) that depends on some external control parameter u (such 
as temperature) that can be continuously varied. In critical 
transitions, @ varies continuously at u, (where it takes a zero 
value) but the derivatives of @ are discontinuous at criticality. For 
the so-called first-order transitions (such as the water-ice phase 
change) there is a discontinuous jump in ¢ at the critical point. 

Although it might seem very difficult to design a microscopic 
model able to provide insight into how phase transitions occur, it 
turns out that great insight has been achieved by using extremely 
simplified models of reality. In this section we introduce the 
most popular model of a phase transition: the /sing model (Brush 
1967; Stanley 1975; Montroll 1981; Bruce and Wallace 1989; 
Goldenfeld 1992; Binney et al. 1993; Christensen and Moloney 
2005). Although this is a model of a physical system, it has been 
used in other contexts, including those of membrane receptors 
(Duke and Bray 1999), financial markets (Bornhodlt and Wagner 
2002; Sieczka and Holyst 2008), ecology (Katori et al. 1998; 
Schlicht and Iwasa 2004), and social systems (Stauffer 2008). 
Early proposed as a simple model of critical behavior, it was soon 
realized that it provides a powerful framework for understanding 
different phase transitions using a small amount of fundamental 
features. 
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The Ising model can be easily implemented using a computer 
simulation in two dimensions. We start from a square lattice 
involving L x L sites. Each site is occupied by a spin having 
just two possible states: —1 (down) and +1 (up). These states 
can be understood as microscopic magnets (iron atoms) pointing 
either north or south. The total magnetization M(T) for a given 
temperature T is simply the sum M(T) = (1/N) Pph S; where 
N = L?. Iron atoms have a natural tendency to align with their 
neighboring atoms in the same direction. If a “down” atom is 
surrounded by “up” neighbors, it will tend to adopt the same “up” 
state. The final state would be a lattice with only “up” or “down” 
units. This defines the ordered phase, where the magnetization 
either takes the value M=1 or M = —1. The system tries to 
minimize the energy, the so-called Hamiltonian: 


1 
Uen (1.3) 
4] 


where J > 0 is a coupling constant (the strength of the local 
interactions) and >? ,; j indicates sum over nearest neighbors.“ 

The model is based on the following observation. If we heat 
a piece of iron (a so-called ferromagnet) to high temperature, 
then no magnetic attraction is observed. This is due to the fact 
that thermal perturbations disrupt atomic interactions by flipping 
single atoms irrespective of the state of their neighbors. If the 
applied temperature is high enough, then the atoms will acquire 
random configurations, and the global magnetization will be zero. 
This defines the so-called disordered phase. The problem thus 
involves a conflict between two tendencies: the first toward order, 
associated to the coupling between nearest atoms, and the second 
toward disorder, due to external noise. 


4 This model can be generalized by including an external field 4; acting on 
each spin and favoring the alignment in a preferential direction. This can be 
done by adding a term — J`; 4; S; to the Hamiltonian. 
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We can see that pairs of units with the same value will 
contribute by /owering the energy whereas pairs with opposite 
values will increase it. Thus changes resulting from spin-spin 
interactions will occur in the direction of reducing the energy of 
the system by aligning spins in the same direction. 

The reader can check that the minimal energy state (the 
ground state) is either all spins up or down (in both cases, we 
have S;5; = 1 for all pairs). However, as soon as we consider 
thermal fluctuations, the system will be able to explore different 
configurations. We will be interested in those more likely to 
occur. In this context, it can be shown that the probability 
P({S;}) of observing a given set {5;} of spins (a state) is: 


oe) 


AT (1.4) 


1 
PS) = oP (- 


where Z is the normalization Z = ts} exp(—E({S;}/&T) also 
known as the partition function.’ Here exp(—E({5;})/&T) is the 
so-called Boltzmann factor. As we can see, increasing T has a 
strong effect: the exponential decay becomes slower and for high 
T the factor approaches 1 for all energies. The flatness implies that 
all states are equally likely to be observed. A typical configuration 
will be a mixed set of randomly distributed spins. 


1.5 Monte Carlo Simulation 


The model runs start with a completely disordered set, where each 
point in the lattice can take up or down states with the same 
probability. Using a fixed temperature T and coupling constant 
J, the new state is obtained by means of a simple set of rules to be 


> The partition function captures all the relevant information required in order 
to recover all the thermodynamic functions, including free energy or entropy, 
among others. 
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applied iteratively using a so-called Monte Carlo method (Landau 
and Binder 2000): 


1. Choose a random spin S;. 

2. Compute the energy change AE associated to the flip 
S; = —S;. 

3. Generate a random number 0 < £ < 1 with a uniform 
distribution. 

4. If E <exp(—AE/&T) then flip the spin.® If not, leave 
it in its previous state. 


5. Repeat [1]. 


The energy difference AF weights the likelihood of the 
transition $; > —S; of taking place. The probability for such 


transition, will be 


W(S 5 ae 1 

( i? )= ep (-25) (1.5) 
if AE <0 and W(S; > —S;) = 1 otherwise. As defined, we can 
see that the temperature is randomizing the system’s behavior (the 
higher 7, the more close to a coin toss we get). This stochastic 
set of rules is known as the Metropolis algorithm (Landau and 
Binder 2000). By iterating these rules, we decrease the energy of 
the system and thus approach the most likely state. In figure 1.5 
we plot the result of computer simulations using a 100 x 100 
lattice. 

We should not be surprised if this caricature of reality fails to 
reproduce anything quantitatively relevant. But let us simulate it 
in two dimensions using the temperature as a tuning parameter. 
Each time step, the magnetization M(T, t) is calculated. Starting 
from a random initial condition, we can compute the order 
parameter, defined here as the average (M(T)) for different 


6 The & parameter is the so-called Boltzman constant, which for convenience 
we take as unity in our simulations. 
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Figure 1.5. Dynamics and spatial patterns in the two-dimensional 
Ising model. Here (a) shows that the transition takes place at 7; = 
2.27. Here a 100 x 100 system has been used. Starting from a 
random initial condition, for each temperature the system is run 
over 5 x 10” steps and afterward M is averaged over 10° steps. Two 
different sets of initial conditions have been used, starting with more 
(filled triangles) and less (open triangles) than half-up spins at time 
zero. In (b) we display a time series of M for T.. 


temperatures. Specifically, the model is run for a large number of 
steps and the average is computed over the last T steps (here we use 
tT=5 x 104, see figure 1.5). Afterward M(7)) is plotted against 
the control parameter 7. In other words, we plot the average: 
1 +t 
(M(T,t)); =- > MAT, t) (1.6) 
Two sets of initial conditions have been used to generate the two 
curves shown here. Open and filled triangles indicate an initially 
higher and lower number of up spins at time zero, respectively. 
A sharp transition takes place at a critical temperature 7}. For 
low and high temperatures, we find what we expected (inset 
figure 1.5a): a homogeneous (ordered system) with |M| ~ 1 anda 
random (M ~ 0) system, respectively. But at criticality the system 
displays wide fluctuations in M, with an average value of zero. 
This is shown by the time series of M(T, t) displayed in 
figure 1.5b. As we can clearly appreciate, wide fluctuations are 
present, far from what would be expected in a system at (or 
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Figure 1.6. Clusters in the Ising model. A cluster involving nearest 
spins having the same orientation (say “up”) is defined by a set of 
spins that are nearest neighbors. Three examples are shown in (a) for 
clusters of size s = 1, s = 2 and s = 4. Here we define the members 
of a cluster by considering only four nearest positions. The cluster 
size distribution at criticality (i.e., for T = T}) is displayed in (b) in 
a log-log plot where data are binned using powers of two. The slope 
T is close to two. 


close) to equilibrium. These fluctuations are a characteristic of 
critical points, but we will ignore them in our treatment of phase 
transitions. Similarly, we could also measure the distribution of 
“islands” of different sizes, composed by groups of spins with the 
same orientation (figure 1.6a). If such distribution is computed, 
it is found that most clusters are composed of only one element, 
followed (with much smaller numbers) by two-spin clusters, and 
so forth. The shape of this distribution is a power law. Specifically, 
if P(s) indicates the frequancy of clusters of size s , this probability 


decays as: 
P(s)~s-* (1.7) 
In a log-log plot, this will approach a straight line (figure 1.6b). 


Power laws are associated to many different complex systems in 
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nature and society (Solé and Goodwin 2001; Ball 2008) and 
indicate the presence of large fluctuations. 


1.6 Scaling and Universality 


We have used the Ising model as an illustration of how the conflict 
between order and disorder leads to long-range correlations in 
equilibrium systems, where a well-defined stationary distribution 
of energy states exists. Other models of equilibrium phenomena 
display similar behavior. In spite of the presence of a characteristic 
scale of local interactions, structures at all scales emerge. But 
the interest of the model resides in its universal character and 
illustrates the astonishing power of model simplicity (Chaikin 
and Lubensky 1995; Stanley et al. 2000). One can determine a 
number of different quantities from the model, such as fractal 
dimensions or correlation lengths. Near the critical point, the 
system displays scale invariance: If we observe such system at a 
spatial scale r, it looks the same when a larger or smaller scale is 
used.” The behavior of these quantities near critical points can be 
fully characterized by a set of numbers, so-called critical exponents 
(Stanley 1972; Goldenfeld 1992). For example, if we measure the 
magnetization M close to criticality, it is found that 


M~(T-T,)? (1.8) 
where B(d = 2) = 1/8 and B(d = 3) = 0.325 in two and three 


dimensions, respectively, to be compared with the experimental 
result for ferromagnetic materials, 8. © 0.316 — 0.327. Similarly, 
the correlation length & scales as 


of = T” (1.9) 


7 Formally, the same statistical patterns are observed when a transformation 
r — Ar is performed (for A > 0). 
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where now v(d = 2) = 1 and v(d = 3) = 0.630, to be compared 
to the real value v, ~ 0.625. The correlation length measures the 
characteristic distance over which the behavior of one element 
is correlated with (influenced by) the behavior of another. Its 
divergence at 7; implies that two distant points influence each 
other and is associated to the emergence of very large clusters, as 
shown in figure 1.6. 

These exponents have been shown to be basically independent 
on the special traits of both interactions and components.® The 
exponents predicted by the model are exact: they fit experimental 
results almost perfectly. Actually, other very different physical 
systems (such as fluids) displaying critical phase transitions have 
been shown to behave exactly as ferromagnets (Yeomans 1992). 
The irrelevance of the details becomes obvious when we check 
that the exact lattice geometry or taking into account interactions 
beyond nearest neighbors does not modify the system behavior. 
Although the value of T} is system dependent, the values of the 
critical exponents appear to be rather insensitive to microscopic 
details, thus indicating that the properties displayed by critical 
points are rather universal. This allows the association of each type 
of transition with a different universality class. The number of such 
classes is limited and thus in principle we can classify different 
phenomena as members of different classes of universality. 


1.7 Mean Field Ising Model 


The set of rules defining the Ising model are simple, so much 
so that it might seem an easy task to predict how the system 
will behave. This is far from true. As discussed at the beginning 
of this chapter, complexity arising from simple interacting units 


8 The concept of universality was originally developed by Leo Kadanoff and 
establishes that near critical points the details of what happens at the small scale 
are largely irrelevant. 
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Figure 1.7. Mean field approach to the Ising model. In (a) a small 
lattice of spins is shown, showing both up and down states. Although 
each spin interacts only with a reduced number of neighbors (eight 
in this case) the mean field approximation (b) reduces to a much 
simpler scenario. A given, arbitrary spin S; would be affected by an 
average field described in terms of the global average spins in the 
system. 


is not simple nor reducible. Solving the two-dimensional Ising 
model was an extraordinary task, finally achieved by Lars Onsager 
in 1944. Among other things he predicted an exact critical 
temperature 7; = 2/log(1 + V2) © 2.27 (Onsager 1944). The 
three-dimensional model was never solved exactly. It has been 
suggested that no such solution actually exists (Istrail 2000). 

In this book we consider the simplest approximation that 
can be made toward understanding phase transitions: the so- 
called mean field theory, which is based on a number of strong 
assumptions and ignores the local nature of interactions. We 
know that each element in the Ising model (or in a real solid) 
interacts with just a few neighbors. We will ignore this and 
assume that each unit “feels” all the others as some kind of global 
field (figure 1.7). This field is nothing but the average value of 
all other spins. In this way, we ignore the spatial and temporal 
fluctuations that we know are present in the real counterpart. This 
is certainly a very strong assumption and as a result it often fails 
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to reproduce some essential quantitative features. Nevertheless it 
captures the /ogic of the system and allows one to predict whether 
phase changes are likely to occur. 

Let us illustrate the approach by solving the Ising model 
(figure 1.6). Consider an arbitrary spin in the system, §S;. 
The average value of this spin can be obtained from (5S;) = 
cer ee S; P(S;) where P(S;) indicates the probability that the 
spin has a +1 or —1 state. Using the Boltzmann distribution, 


we have 
al [2 Çs >j s)| (1.10) 
2 s;=—1,41 XP [-2 (JS: Ši s;)| 
and the average spin value now reads: 
exp |- @ Dag s;)| -P E (Uy dj s;)| (1.11) 


a ES +0 [- 75,5) 


the previous expression can be transformed using the equivalence? 
tanh(ux) = [exp(— ux) — exp(ux)]/lexp( ~ux) + exp(ux)] 
which for our case gives: 


(S;) = —tanh |B [JX S; (1.12) 
j 


This function is known to be bounded by two assymptotic limits, 
namely tanh(ux) > +1 for x — +00. 


Let us look at the following sum: H; = o2. J5;) S; which 
is nothing but the (local) energy of interaction associated to the 
S; spin. Now let us make the mean field assumption: replace the 


? The hyperbolic tangent tanh(x) and other so-called hyperbolic functions 
are the analogs of the standard trigonometric functions, with several desirable 
properties, and widely used in many areas of physics and mathematics. 
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previous spins by their average values, namely (S;) = yy S;/N, 
meaning that the exact set of states surrounding S; is replaced by 
the average over the lattice. Now we have 


= 115) > 7 |S (1.13) 
J 


If the number of neighbors is g (four in two dimensions, eight in 
three, etc.) we can write X` j J =J¢ and the average value of the 
spin is now: 


(S;) = tanh [B (4J (S;))] (1.14) 


Since we defined the magnetization M as the average over spins, 
and using B =1/kT we have a mean field equation for the 
magnetization in the Ising model: 


= 4r 
M=tanh | (2) M (1.15) 


where we used M = (s;) = (s ;). For convenience (as we will see 
below) we will use the notation 7; =qJ//k. In figure 1.8a we 
plot the function yı = tanh(M) for different values of u = T /T 
which acts as an inverse of the temperature. We also plot the linear 
function y2 = M. The potential equilibrium states of our system 
would correspond to the intersections between these two curves, 
that is, yj = y2. As u is increased (and thus the temperature 
reduced) the steepness of the nonlinear function y; increases close 
to M=0. At high T (low u) only one intersection is found, 
corresponding to the zero magnetization state. Once we have 
T < T. three intersections occur. 

If we are interested in studying the behavior of this system close 
to the phase transition, where M is (on average) small, we can use 
the expansion tanh x © x — x?/3 +--- By ignoring all terms of 
order larger than three, we obtain an approximate expression for 
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Figure 1.8. The mean field theory for the Ising model correctly 
predicts a phase transition close to a given critical temperature T}. 
In (a) we show the shape exhibited by the function y = tanh(u M) 
against the magnetization M for different u values. The line y = M 
is also shown, together with the intersections with the curves, 
indicating solutions to the equation M = tanh(u M). In (b) the 
resulting transition curves are shown. 


the magnetization: 


T: DE 
M=| Z] M=-| ZM (1.16) 
T SL 
Which values of M(T) are compatible with this equation? One 
solution is the trivial, zero magnetization state, that is, M(T) = 0, 
and the other possible values follow 1 = (7, /T)—(TZ./T)°/3M. 
Two symmetric solutions are obtained, namely: 


T ‘die 
M,(T) =+y/3 (=) h = J (1.17) 


which exist only if T < T.. These solutions are shown in figure 
1.7b using the reduced temperature T/ 7; as our control param- 
eter. As we can see, the mean field model properly captures the 
key behavior of the magnetic phase transition. Several estimations 
of relevant thermodynamic quantities can be derived from this 
approximation (see Christensen and Moloney 2005). 
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Since the mean field approach ignores spatial and temporal 
correlations (and fluctuations) it must somewhat fail to be a 
correct description of real systems. As a consequence, predicted 
scaling laws and other quantitative measurements will depart 
from the observed.'° However, this theoretical approach gives in 
general the correct qualitative predictions of the phase diagrams, 
particularly as we increase the range of the interactions. We 
can imagine this as follows: in a d-dimensional square lattice, 
the number of nearest neighbors is 2d, which means that, as 
d grows (i.e., as more neighbors are included) each element 
exchanges information with a larger fraction of the system. As 
a consequence, the effects of spatial degrees of freedom and 
fluctuations become less and less important. In fact, statistical 
physics analyses have shown that there is a critical dimension 
d, above which fluctuations become unimportant and the MF 
approach becomes a correct description (Goldenfeld 1992). 


In summary, as the number of potential interactions grows, mean 
field predictions become more and more exact. Moreover, as 
pointed out by Yeomans (1992) the mean field approximation 
provides a feasible approach to complicated problems and, in 
many cases, is the only one. Moreover, it often provides the first 
step toward more accurate models. For all these reasons, this is 


the path taken here. 


1.8 Nonequilibrium Transitions 


Is the Ising model our model of reference for complex systems? 
I think it provides an elegant example of how global complexity 


10 The previous formalism can be generalized by using a continuous approx- 
imation where the local states are replaced by a field @(x) now x indicating 
the coordinates in a continuum space. The so-called field theory approximation 
allows us to explore the effects of spatial degrees of freedom and is the 
cornerstone of the modern theory of critical phenomena (Kardar 2007). 
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emerges out from local interactions. It also illustrates the nontriv- 
ial phenomena arising close to critical points. But the hard truth 
is that most systems in nature and society are far from equilibrium 
and no energy functions are at work (Hinrichsen 2006). In 
particular, equilibrium systems obey the so-called detailed balance 
condition. Specifically, if the system is described in terms of a 
set of possible states, A= {a} such that a € A occurs with a 
probability Py, for each pair a, a’ € A, a transition probability 
æla — a’) can be defined. The detailed balance condition is 


defined as follows: 
P ola > a’) = Pyol(a’ > a) (1.18) 


which can be easily interpreted: the probability flows from one 
state to another, in both directions, and cancels the other out. 

Nonequilibrium systems are not determined (solely) by exter- 
nal constraints, and violate the detailed balance condition. As a 
given control parameter is changed, a given stationary state can 
become unstable and be replaced by another. Since detailed bal- 
ance is not operating, their macroscopic properties are no longer 
obtainable from a stationary probability distribution. Besides, in 
systems out of equilibrium, no free energy can be defined. In 
equilibrium systems, the phase diagram can be determined on the 
basis of the free energy but this approach becomes useless when 
we move far from equilibrium. 

The difference is particularly dramatic when dealing with 
systems involving absorbing states (Hinrichsen 2000; Marro and 
Dickman 1999). Absorbing configurations can be reached by the 
dynamics but cannot be left. As an example, a population can go 
extinct once a given threshold is reached, but there is no way to 
get out from extinction. Not surprisingly, extinction and collapse 
will be two important examples of nonequilibrium transitions 


explored in this book. 


2 


STABILITY AND INSTABILITY 


2.1 Dynamical Systems 


Our first basic approach to the analysis of phase changes in 
complex systems involves the consideration of a fundamental 
property: stability. Our perception of the stable and the unstable 
is filtered by our perception of change. But in general terms 
it can be said that stable things are those that do not vary 
over long time scales, whereas instability is tied to change. 
Complex systems maintain a stable organization thanks to a 
balance between internal order (which needs to be preserved) and 
flexibility (adaptability). 

A mechanical analog of the basic ideas to be discussed below 
is displayed in figure 2.1. Here we have a surface on top of which 
one can place marbles. Although a marble could be left on top of 
a peak and remain there, any small perturbation would force the 
ball to move away (and thus the peaks are unstable). Instead, the 
bottom of a valley is a stable point, since small perturbations of 
the marble only change its position for a while, until it goes back 
to the bottom again (to the stable equilibrium state). Marginal 
points are represented by locations on flat domains, where a small 
movement to a neighbor location just leaves the ball in the same 
place, without moving away or back to the initial position. 
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Unstable 


Marginal 


Figure 2.1. The three types of (linear) stability, through a simple 
mechanical analogy. Here three marbles are at equilibrium at three 
different points in the landscape (From Solé and Bascompte 2006). 


The simplest dynamical system that we can consider is de- 
scribed by a /inear differential equation. Imagine a population of 
bacteria replicating at a given rate. Let us indicate as u the rate of 
replication of individual bacteria. Clearly, the more bacteria we 
have in our Petri dish, the faster the production of new bacteria. 
If we indicate by means of x(t) the number of bacteria at a given 
time ¢, after a small time At we would have: 


x(t + At) = x(t) + ux(t)At (2.1) 


This discrete equation just tells us that the process is incremental 
and proportional to the existing population. The last term on the 
right-hand side indicates that the longer the interval Az, the larger 
the number of newborn bacteria. 

Although this equation is useful, it is rather limited: we would 
like to know how many bacteria are expected to be found after a 
given time, assuming that the initial number x (0) at a given zero 
time is known. In order to find this, we need to transform the 
previous difference equation. 

The time interval Aż is an arbitrary one, and we can make it 
as small as we want. By rearranging the previous linear equation, 
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we have: 
x(¢ + At) — x(t) _ 
At E 


Using x = x(t), and considering the limit when At — 0, it can 


px (t) (2.2) 


be transformed into a differential equation: 
— = ux (2.3) 


where dx /dt is the derivative of x against ż, measuring the speed 
at which the population grows with time. As we can see, this speed 
is proportional to the current population, increasing linearly with 
x. For this reason, we define this model as a one-dimensional linear 
dynamical system. 


2.2 Attractors 


The linear differential equation can be used to find the explicit 
form of the x(ż) time dependent population dynamics. Solving 
the differential equation in this case is an easy task. We just 
separate the variables, that is: 


a = udt (2.4) 
x 


and perform the integration in both sides. Assuming that at time 
t = 0 the initial population was x9 and that the population size 
at any arbitrary time ¢ > 0 is indicated as x(t) we have: 


x(t) t 
f k = f udt (2.5) 
x0 x 0 


which gives the solution 
x(t) = xoe* (2.6) 


The predicted behavior indicates that an explosive (exponential) 
growth will take place, with rapidly increasing numbers of indi- 
viduals as time proceeds. 
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A very different scenario is provided by the linear decay defined 
by the differential equation: 

dx _ 

P 


where the speed of change is now a negative linear function 


—ux (2.7) 


of population size. This equation is also common in different 
contexts, such as radioactive decay:! if x indicates the number of 
atoms of a given radioactive material, they disintegrate at a given 
rate u. The corresponding solution is now 


x(t) = xoe (2.8) 


where xo would represent the total number of unstable atoms at 
the beginning of the experiment. In this case, when we consider 
the limit at infinite time, we obtain 


lim x(t) = lim xe = 0 (2.9) 


t>co t—> 00 


since asymptotically all atoms will eventually decay. 

Sometimes, we can reduce the complexity of a given problem 
involving several variables to a single-equation model. Let us 
consider an example that has some interesting features. The 
problem is illustrated in figure 2.2, where we display a schematic 
diagram of two connected compartments where a gas is confined. 
Initially, all molecules are placed in the left compartment. The 
two boxes are connected through a valve of fixed diameter that is 
open at time zero. Once open, the gas particles can flow between 
the two compartments. What is the dynamics of gas exchange? 
We know that eventually the molecules will be equally distributed 
between both compartments and such a result can be easily shown 
by means of an appropriate use of mass conservation. 


! Radioactive decay is a stochastic process: any given radioactive atom of a 
given isotope can disintegrate at any time with a given probability. However, 
taken together and considering a large number of atoms, the statistical behavior 
of the material is well described by a deterministic linear model. 


Figure 2.2. An example of a simple physical system where conser- 
vation laws permit reduction of its mathematical description. Here a 
two-compartment system is considered, with a gas initially located in 
the left compartment (a) that diffuses through time until it reaches 
equilibrium (b). The system is closed and the total number of 
particles is thus conserved. 


Let us call C the initial concentration of molecules in the first 
compartment. Once the gas starts flowing after opening the valve, 
the concentrations in each box (indicated as x and y) will change 
in time, according to a system of two differential equations: 


dx 

T = p(y — x) (2.10) 
T 2.11 
p oN (2.11) 


which are clearly linear.? The parameter p is a so-called diffusion 
constant and weights how easily particles move from one com- 
partment to the other. 

Since we assume a closed system, we have x + y = C and thus 
we have: 


dx 
g MC — 2x) (2.12) 


which gives an evolution equation: 


x(t) = a =g t] (2.13) 


? In general, a system of 7 linear differential equations will read as dx;/dt = 
j UGX) ford, fl ccc, n and being œ;; a matrix of real coefficients. 
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As expected from our intuition, the final state is determined 
by the limit: 


C C 
lm [=e t] = — (2.14) 


tooo 2 


N 


and thus the equilibrium state is x* = y* = C/2: the gas uni- 
formly distributes over space. Any initial condition will eventually 
end up in this equilibrium state. The example thus illustrates one 
possible scenario where dimensional reduction can be obtained by 
means of conservation. 

We must mention here that linear differential equations are 
seldom a good description of complex systems, but as will be 
shown in the next section, the previous examples will be useful 
in developing a general approximation to nonlinear dynamics. 
Moreover, although we have solved analytically the previous 
models, finding an explicit form for x(t) is in most cases impos- 
sible. However, since we are mainly interested in the observable 
stationary states, we will use a different strategy of analysis. This 
strategy makes use of a general technique that determines the type 


of stability displayed by each fixed point. 


2.3 Nonlinear Dynamics 


In general, we will consider dynamical systems described by 
means of nonlinear differential equations. The general expression 
in our one-dimensional representation would be 


d 
on = f(x) (2.15) 


where f;,(x) is some type of nonlinear function of x (such as 
x*, x3, cos(ux) or 1/(u +x)). An alternative notation commonly 


used is: 


x = file) (2.16) 
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where x indicates derivative of x over time. The function is 
assumed to be continuous and with continuous derivative. The 
u symbol is used to remind us that in general the dynamics will 
depend on one or several parameters. 

As an example of a well-known nonlinear system, consider the 
logistic growth model (Pielou 1977; Case 1999). If x indicates 
the population size of a given species growing on a limited set 
of resources, its time evolution can be described by means of the 
following deterministic equation: 

d 


x x 
Se = ful) = x (1 = =) (2.17) 


where u and K are the growth rate and carrying capacity of the 
population, respectively. We can easily see that when x is very 
small compared to K it can be approximated by dx/dt © ux, 
which means that an exponential growth should be expected. 
However, the growth will be limited by available resources, and 
for x ~ K (close to the limit) the last term in the right-hand side 
of the logistic equation will go to zero, and thus the speed of 
population growth will decrease to zero too. These patterns can 
be studied by solving the differential equation analytically.? If 
x = xo at t = 0, it can be shown that 


K 
a 1+ [l enht 


Xo 


(2.18) 


and thus x(t) > x* = K as t— oo. In other words, any initial 
condition xp >0 eventually converges to the carrying capacity 
x* = K, consistent with our intuition. In figure 2.3a we show sev- 
eral examples of the trajectories followed by the logistic dynamics, 
as described in equation (17) using u = 1 and K = 2 and 


> Once again, we can separate variables and solve the integrals f dx/(x(1 — 
x/K)) = f udt. This is done by decomposing the x-dependent fraction as a 
sum 1/(x(1 — x/K)) = A/x + B/(1 — x/K) and then determining the values 
of A and B consistently and integrating separately the two terms. 
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Figure 2.3. Solutions of the logistic differential equation, for u = 
1.0 and K = 2. In (a) several initial conditions are used, all of them 
converging to the steady state x* = K, which acts as an attractor. In 
(b) we show a single solution, with x(0) = 0.01 in linear-log scale. 
The dashed line indicates the exponential law x(t) © x(0)e"* which 
is valid for short times. 


different initial conditions. We can see how all of them approach 
the carrying capacity as time proceeds, thus indicating that x = K 
is a stable state. In figure 2.3b we show one of these trajectories 
(starting from x(0) = 0.01) in linear-log scale. In such a plot, 
exponential solutions appear as straight lines. As we can see, 
this is the case for the first part of the dynamics, as expected: 
The exponential solution obtained from the linear approximation 
dx/dt ~% ux is also observable (dashed line) showing that the 
linear approximation is good for small time scales. Since all 
trajectories starting from x9 > 0 become eventually “attracted” by 
x*, we call it an attractor of the dynamics. This is also the case 
for other models, in which stable points describe the long-term 
behavior of the solutions. 

Instead of finding the analytic solution of the differential 
equation (which is not possible most of the time) we can study the 
behavior of the solutions close to the fixed points. More precisely, 
we want to determine how to classify different types of fixed 
points in terms of their stability. In order to do so, we will study 
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the behavior of the system close to the equilibrium points, using a 
linear approximation to the nonlinear dynamics. This defines the 
so-called linear stability analysis (LSA). 


2.4 Linear Stability Analysis 


The first step in LSA is to identify the set 7, of equilibrium (fixed) 
points. Since equilibrium will occur whenever dx/dt = 0 (no 
time changes take place) this is equivalent to finding those x* 
such that the function f;,(x*) is zero. The resulting set of points 


will be: 
Ty = {x* | fale) = 0} (2.19) 


From our definition, if x(0) = x*, it is true that x(0o) = x*. 
But obvious differences appear when comparing the behavior of 
different fixed points. The basic idea is illustrated using again 
the mechanical analog with marbles shown in figure 2.1. The 
three of them are at equilibrium. However, a small displacement 
from their equilibrium state will affect them in very different 
ways. The stable equilibrium point will be recovered after some 
transient time (as it occurs with a pendulum). The unstable one 
will never recover its original position: the initial perturbation will 
be amplified and the marble will roll away. A third, apparently less 
interesting possibility is the so-called marginal state: on a totally 
flat surface, a marble displaced to another close position will just 
stay there. 

The nature of the different points x* € 7r, is easily established 
by considering how a small perturbation will evolve in time 
(Strogatz 1998; Kaplan and Glass 1999). Consider x = x* + y, 
where y is a small change (such as adding a small amount of 
individuals to the equilibrium population). Using the general 
form of the one-dimensional dynamics given by eqn (15) and 
since x* is a constant, the dynamics of the perturbation y is 
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obtained from: 


dy 
= falx* + y) (2.20) 
Assuming that y is very small, we can Taylor expand‘ the previous 
function: 
dy 


oh A y (2.21) 


Z = f) + |e trl de 


If y is very small, we can safely ignore the higher-order terms (with 
y’, y? .). Since fu (x*) = 0 (by definition) and provided that its 


first derivative is non-zero, we can use the linear approximation: 


2 = Any (2.22) 
where À, is defined as: 
àu = K (2.23) 
or alternatively as 
Àu = E) (2.24) 


The linear system is exactly solvable, with an exponential 
solution: 


y(t) = y (0) (2.25) 


This solution immediately provides a criterion for stability: if 
hy < 0, then y — 0 and thus the system returns to its stable 
equilibrium. If à, > 0, perturbations grow and thus the point 


4 In general, a Taylor expansion of a given function about a point x* is given 
by an infinite series, namely f(x) = f(x*) + f'(x*)(x — x*)/L 4 f(x") (x — 
x*)?/2!4+---.In our notation, we have y = x — x* which gives the size of the 
perturbation relative to x*. 
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Figure 2.4. The fixed points of a given dynamical system can be 
easily obtained by plotting the function f (x) against x and finding 
its intersections with the horizontal axis. Here the logistic curve is 
shown for u = 1 and K = 5. Stable and unstable fixed points are 
indicated as filled and open circles, respectively. 


is unstable. Both possibilities are separated by a marginal state 
defined from àf, = 0. Thus just looking at the sign of A, we 
can immediately identify the presence of stable and unstable 
states. Using the condition à, = 0 we will also characterize the 
boundaries separating stable and unstable patterns of behavior. 

This finding can be geometrically described by plotting the 
function f,,(x) against x (figure 2.4). The intersections of f with 
the horizontal axis correspond to the fixed points, and the slope 
of the curve near them is related to their stability. At unstable 
points, the derivative of f,(x) is positive, which means positive 
slopes. Instead, stable points will have negative slopes. By plotting 
this picture, we can immediately identify how many fixed points 
are present and their stability. 
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We can illustrate these general results by means of the logistic 
model just described. Here we have a set of fixed points defined 
by 2, = {0, K} and the stability of each equilibrium point is 
determined by the sign of 


p) x 
hy =U (: = - ) (2.26) 


In this case we have we have 4(0) = u which is always positive 
(since the growth rate is a positive parameter) and thus x* = 0 
is unstable. For the second point we have A(K) = —2 which 
is always negative and K is consequently always stable, consistent 
with our intuition. 

Stability theory, as defined here, will be our first tool in 
the study of mean field models of phase changes. The second 
ingredient has to do with stability changes: how and when a 
given stable state becomes unstable, being replaced by a new, 
qualitatively different type of state. 


3 


BIFURCATIONS AND CATASTROPHES 


3.1 Qualitative Changes 


As discussed in chapter 1, our interest here is understanding the 
presence of different types of behavior defining different phases. 
As it occurs with physical systems, phases can be defined and 
their boundaries characterized. The stability analysis presented in 
the previous chapter allows exploring changes in the behavior of 
complex systems in terms of changes in their stability. 

We have already mentioned a few examples of phase changes in 
the first chapter. Here our goal is to provide a simple but useful 
mathematical approach to characterize the problem. Following 
the formalism introduced in the previous chapter, we will con- 
sider models that allow describing a given phase change in terms 
of one-dimensional differential equations. Such a simplification, 
as already discussed, implies a number of limitations. But it also 
helps understanding how we can capture some key aspects of 
complexity by properly choosing the right observable. 


The key concept presented here is the idea of bifurcation (Jackson 
1990; Nicolis 1995; Kuztnesov 1995; Strogatz 1998), a term 
attributed to Henri Poincaré. A bifurcation is a qualitative 
change in the dynamics of a given system that takes place under 
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continuous variation (tuning) of a given parameter. Bifurcations 
will be the analog to phase transitions under the continuous, 
dynamical systems view taken here. Mathematically a bifurcation 
implies the emergence (and disappearance) of new solutions. 
In the real world, it relates to a qualitative change when, for 
example, information, viruses, or fires propagate (what thresholds 
need to be overcome) or a shift occurs in the structure of an 
ecosystem. An example of the latter would be the transition 
from vegetated to desert habitats. Moreover, some specially 
relevant bifurcations (see next section) imply the appearance 
of alternate states. Such bifurcations introduce an interesting 
element: the possibility of path-dependent phenomena. In these 
cases, relevant to both biological and technological evolution 
(Arthur 1994; Gould 1992), the initial condition plays a key 
role. Small differences in the initial state can actually lead to very 
different outcomes, thus introducing a role for history. 


3.2 Symmetry Breaking 


In many situations, a given system is forced to make a choice 
between two potential candidate states. These two states can be 
equally likely chosen, and the final choice is often decided by 
means of some historical contingency. Once again, the situation is 
illustrated by means of a mechanical analogy, shown in figure 3.1. 
In this case, a single marble is located on top of a ridge separating 
two identical valleys. Any small change in the marble’s position 
will force it to roll down toward the nearest stable state (the 
bottom of the nearest valley). The direction of such initial 
perturbation plays a relevant role here, since only one of the two 
possible choices can be made. The initial symmetry is broken. 
Such symmetry breaking (SB) mechanism allows us to create 
structures and pervades the generation of patterns in many 
situations in biology, chemistry, and physics (Haken 1975, 1996; 
Nicolis and Prigogine 1989; Strocchi 2005). The two alternative 
orientations of spins displayed by the Ising model (figure 1.5) is a 
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Figure 3.1. Symmetry breaking: a marble can roll down on a 
landscape where two alternative (stable) valleys (S) can be reached. 
The valleys are identical and located symmetrically from the original 
position of the marble, which is in an unstable state (U). After Solé 
and Bascompte, 2006. 


good illustration of this phenomenon. By depending on the initial 
balance between up and down spins, the magnet evolves toward 
either a positive or a negative global magnetization, respectively. 

The SB phenomenon is very well illustrated by the following 
dynamical model: 


d 

E = px — x’ (3.1) 
Which appears in many different contexts. The set of fixed points 
of this system is easily determined. Together with the trivial fixed 


point x* = 0 we have two symmetric points, namely 


= 24/0 (3.2) 
The stability of these fixed points is determined by: 
Ay = u 3(x*)? (3.3) 


For x* = 0, we obtain 4,,(0) = u: This point is thus stable if 
u < 0 and unstable otherwise. For the two other fixed points, 
we have now à (xž) = —2u which means that they both exist 
and are stable for u > 0. These fixed points and their stability 
are easily determined by drawing f;,(x) and looking at its slope 
around the intersections with the x-axis (figure 3.2). 
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Figure 3.2. The presence and type of equilibrium points in the 
dynamical system dx /dt = ux—x?’. For u > 0 (left) three crossings 
occur between ux — x? = 0 and the x axis, thus indicating the 
presence of three fixed points, one of them unstable (open circle) and 
two stable ones (filled circle). For u < 0 (right) only one crossing is 
present, indicating a single (stable) fixed point. 


In figure 3.3a we show the trajectories x(t) for this system, 
using a negative value of the bifurcation parameter u and starting 
from different initial conditions. The corresponding solutions are 
shown in figure 3.3b for a positive value of u. These curves 
have been calculated by numerically integrating the differential 
equation using the so-called Euler’s method.! As predicted by the 
LSA, solutions converge to x* = 0 in the first case while they split 
into two different attractors for the second case. 

We can integrate all the previous results by using a simple 
diagram. Such a bifurcation diagram provides a picture of the 
available equilibrium states and their stability as the bifurcation 
parameter ju is tuned. In figure 3.3c the bifurcation diagram for 


! Euler’s method is the simplest algorithm that allows solving of a differential 
equation. It is based on the discrete approximation to the differential equation 
already discussed in the introduction section. Given an initial condition x (0) and 
fixing the time interval Aż (typically a small value) we compute x(Az) by using 
x (At) = x(0)+ ff. (0) Az, and the next value from x (2At) = x(At)+ fi. (At) At 
and so on. 
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Equilibrium points (x) 


Figure 3.3. Symmetry breaking in a one-dimensional dynamical 
system, as defined in equation (1). In (a) we show several numerical 
solutions for this model using u = —0.5 starting from different 
initial conditions. The same is shown in (b) for u = 0.5. The 
corresponding bifurcation scenario in displayed in (c). 


this model is shown. It is constructed by representing all the fixed 
points x* € z,, for each value of u. Stable fixed points are joined 
by continuous lines, whereas unstable ones are located along 
discontinuous lines. We indicate by means of arrows the direction 
of trajectories on the x-space toward stable states. The bifurcation 
diagram also allows to define the boundaries and behavior of the 
two qualitatively differentiated phases displayed by this system. 
It is worth noting the similarities between this example and the 
mean field treatment of the Ising model described in chapter 1. 

The qualitative properties of the SB phenomenon can also 
be captured by a mathematical counterpart of the surface-and- 
marble metaphor. This is achieved by using a so-called potential 
function ®,,(x), defined from the relation 

dx d®,,(x) 


dt dx a 


where ®,,(x) would be computed from 


(x) = = f file == fras (3.5) 
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Figure 3.4. The potential ®,,(x) for the system described by equa- 
tion (1) using: (a) u = —1.5, (b) u = —0.4 and (c) u = 1.5. The 
location of the stable and unstable fixed points is indicated by means 
of empty and filled circles, respectively. For y > 0 two alternative 
states are simultaneously available, symmetrically located on both 
sides of x* = 0. Any small perturbation of the zero state will decide 
the chosen equilibrium among the two, thus breaking the symmetry 
of the system. 


It is usually said that a dynamical system that can be written 
as (3.4) derives from a potential ®,,(x). The existence of a 
potential function will depend on the integrability of f, (x) but 
for smooth enough functions? is typically guaranteed. The logic 
of this description is clear if we notice that the same points that 
give dx /dt = 0 will be those such that d®,,(x)/dx = 0 and thus 
corresponds to the extrema of ®,,(x). Moreover, since we have 
d?®,,(x)/dx? = —df,,(x)/dx = —,, the sign of the second 
derivative will be correlated with the type of equilibrium point. 
Specifically, the minima and maxima of ®,,(x) will correspond to 
stable and unstable points, respectively. 
For the previous example, we have 


xí 2 


O,(x) =F — a 3.6) 


? This essentially means that the function has no singularities, which occurs 
when f and its derivatives are continuous. 
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and its shape is displayed in figure 3.4 for two values of u. When 
the bifurcation parameter is negative, this function has a single 
minimum at the bottom of the valley whereas when positive the 
situation dramatically changes: two valleys are now present. Two 
alternate states are available, and the previous stable point is now 
at the top of a ridge. From there, tiny differences will decide the 
final fate of the system. 

Although this book does not consider models including ex- 
plicitely random perturbations (the Ising model is our exception) 
it is worth mentioning such effects in the context of symme- 
try breaking. Instead of making an analytic derivation, let us 
use a simple mechanical analog to illustrate the point. This is 
summarized in figure 3.5, where we show a marble rolling at 
the bottom of a surface. This would roughly correspond to the 
potential function of our previous model (although an additional 
dimension has been added) below but close to the bifurcation 
point. Close to this point, as shown in figure 3.4b, the bottom 
of the valley flattens and the flatness increases close to criticality. 
Now imagine our marble rolling at the bottom of this lanscape, 
sometimes quickened by small, external perturbations. Since the 


Figure 3.5. Fluctuations close to the symmetry breaking transition. 
Here we represent a close-up of two potential functions ¢,,(x) with 
a shallow (a) and a shallower (b) bottom. The flatter the surface, the 
wider will be the wandering of the marble around the equilibrium 
point. 
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valley is rather flat, any small perturbation can drive the marble 
far away from the equilibrium point. As a consequence, the 
marble might spend long periods without even approaching 
the equilibrium state. The closer we are to the critical point, the 
more dramatic become these fluctuations. This has important 
consequences for the relaxation time required to achieve the 
steady state, as discussed in section 3.4. 


3.3 Catastrophes and Breakpoints 


The previous example illustrating the SB phenomenon is known 
in physics as a second order phase transition. In the context of 
the physical theory of phase transitions, such a phenomenon 
takes place as the control parameter u changes and some order 
parameter (our x) moves continuously from one phase to the 
other. However, in some cases a slow change in the control 
parameter leads to a sudden shift in the system’s state or, as 
it is sometimes dubbed, a catastrophe (Thom 1972; Zeeman 
1977; Poston and Stewart 1978; Deakin 1980; Gilmore 1981; 
Arnold 1984; Jackson 1991). This type of phenomenon occurs 
(as discussed in chapter 1) in the water-ice transition at sea 
level pressure: as temperature (our control parameter) decreases 
below the freezing temperature, water displays a transition to the 
solid state. At the microscopic scale, the liquid water with all its 
molecules dancing around experiences a deep geometrical change, 
being that now the molecules are arranged in an ordered manner. 

Rapid shifts from a given state to a new one, separated by 
a finite gap, are known to occur in many different contexts, 
from ecology (Sheffer et al. 2001; Scheffer and Carpenter 2003) 
to climate (Foley et al. 2003) and social dynamics (Ball 2006). 
The changes can be dramatic and an example is illustrated in 
figure 3.6a, where we return to the example discussed in the 
introduction to this book. What is displayed here is a time series 
of the amount of dust measured in ocean core sediments over 


Terrigenous sediment (%) 
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Figure 3.6. Catastrophic shifts: In (a) we display the time evolution 
of the dust measured from sediments at a given location (redrawn 
from Scheffer et al. 2001). This measure provides a surrogate of 
the area covered by terrigenous (wind-eroded) soils. About fifty- 
five hundred years ago, a sudden drop took place (see text). The 
mechanical analog of such catastrophic shifts is given in (b). 


a long period of time, spanning nine thousand years to the 
present. We can see that the amount of dust decreased slowly 
until a well-defined drop occurred around fifty-five hundred years 
ago. The measurement was obtained near the African coast in 
the northern hemisphere, and the shift is actually associated 
with a large-scale, abrupt change (perhaps occurring over the 
course of a century) from a vegetated to a desert state in the 
geographical area that is today known as the Sahara and Sahel 
deserts. Despite that the shift from vegetation to desert was sharp, 
it was associated with continuous variations in the amount of 
solar radiation (deMenocal et al. 2000). In other words, smooth 
changes in external drivers of vegetation dynamics led to a sudden 
change in the ecosystem organization. This situation might be 
repeating itself currently in semiarid ecosystems (Kefi et al. 2007; 
Scanlon et al. 2007; Solé 2007; Kefi 2008; see also chapters 12 
and 13). A mechanical analog of this phenomenon is shown in 
figure 3.6b. Here a marble is shown on top of a folded surface. 
If it rolls down, it will move smoothly until the fold is reached. 
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Afterward, it will fall into a new area, separated from the previous 
by a well-defined gap. Is there a mathematical counterpart to this 
mechanical analogy? An example will help us see that indeed this 


is the case. 
Consider now a slight generalization of our previous model, 
given by: 
dx 3 
—=atpx—x (3.7) 
dt 


Due to the presence of the new parameter œ, the analysis of the 
system becomes more complicated: the set of fixed points 7T, is 
now obtained from œ + jx* — (x*)? = 0 and the general form 
of them is rather involved. We know, however, that either one 
or three solutions will be possible. Instead of finding analytically 
these points, we will use a shortcut. The fixed points are found at 


the intersection of the curve 
yı = ux — x? (3.8) 


and the constant line yz = —a. As it happened with our previous 
system, we have two very different qualitative behaviors. For 
u < 0 a single fixed point is obtained, whereas for positive 
values (an example is shown in figure 3.7a) there is a domain 
of a values, —a, < œ < +a, where three simultaneous fixed 
points are present. As we move the constant line y2, we can see 
that it crosses the y; curve in different ways. For some values only 
one crossing is allowed (giving a single stable point) whereas for a 
domain of values three crossings occur. The impact of these folded 
lines on the system’s dynamics is illustrated in figure 3.7b, where 
we show that the change in the equilibrium states as @ is linearly 
and continuously tuned (inset c). Starting from œ = 0.5 and the 
corresponding value of x¥, we can see that once œ < a, = 0.25 
a dramatic shift takes place, with the state jumping to the lower 
branch x*. This perfectly illustrates the basic principle that for 
these types of systems a slow change of parameters does not always 
produce a slow dynamical response. 
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Figure 3.7. Catastrophes in a simple nonlinear model. In (a) we plot 
the two curves yı = ux — x? (for u = 0.75) and y2 = —a (dashed 
lines) whose intersections define the possible fixed points. In (b) we 
can see how a continuous change in a control parameter can trigger 
a sudden shift (b). Here œ is varied linearly over time (inset c). 


In order to determine the boundaries of the three-state 
domain, which now depend on two parameters, we need to 
compute @,. First, let us find the value x = x, where the extrema 
of ye = wx, — x? occurs. It is obtained from the condition 
dy/dx = 0 and is given by: 


Xe = +E (3.9) 


and the corresponding y, will be the value @,, namely 
2 
4+4% E 
3 y3 


Using the previous results, we can determine the domain D 


(3.10) 


Ao = 


in the (u, œ) parameter space where three fixed points exist. This 
is now shown in figure 3.8a, where the two possible phases are 
indicated. We also show (figure 3.8b) the corresponding surface 
with the location of the fixed points x* = x*(u, œ), as defined by 
the implicit function a + ux* — (x*)? = 0. We can appreciate 
a remarkable shape in the folded surface, not very different from 
the one we already displayed in our mechanical analog. What is 
the importance of this shape? 
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Figure 3.8. Domains of multistability in dynamical systems. In (a) 
we show the associated domain of the (u, œ)-parameter space for 
the previous example. The cusp shown in gray indicates the domain 
where three solutions are possible. In (b) the corresponding folded 
surface associated with the possible fixed points is also shown (see 
text for details). 


A very important consequence of this folded surface is the 
presence of a memory effect or hysteresis. The surface x*(a, m) 
is shown in figure 3.8b, where we start from a given stationary 
state 1 on the surface and continuously tune @ in such a way 
that the system changes continuously on the surface as the control 
parameter is varied. After reaching point 2 the state jumps into the 
lower part of the surface (at 3) and if we keep reducing œ a new 
state is reached under continuous modifications (4). However, if 
we start now on the lower branch again from 3 and now increase 
a, we do not jump back to 2 but instead move on to the lower 
branch until the fold is reached again, this time at point 5. A new 
jump occurs, but the trajectory followed is totally different. 

The sudden changes associated to these transitions can also be 
seen by looking at the corresponding potential function, which 
now reads: 


2 x4 


(x) = —æx = u% + Fa (3:11) 
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By plotting the potential using different values of the bifurcation 
parameters, we can easily appreciate an important difference with 
respect to the potential functions involved in SB phenomena. 


3.4 Critical Slowing Down 


An important property displayed by systems exhibiting 
symmetry-breaking bifurcations is the presence of a divergence in 
the relaxation time close to criticality. Specifically, as the critical 
point is approached, the characteristic time needed to reach the 
equilibrium state rapidly grows. As we have already mentioned, 
when the system exhibits stochastic fluctuations, the variance of 
these fluctuations close to criticality will increase and actually 
exhibit a marked maximum when looking to experimental sys- 
tems (Haken 1977, 1996; Nicolis and Prigogine 1989). Here we 
derive the time required to approach the stationary state asuming 
deterministic dynamics. 

We can illustrate this phenomenon by using the previous 
model (3.1) and calculating the transient time Zj, required to 
go from a given initial condition x(0) to a final state x(Zj,) as 
a function of the control parameter u. In general, such transient 
time will be obtained from: 


‘a x(Ty) 7 
f j= f Eo de (3.12) 
0 x(0) 


which needs to be split into two terms, separating the two 
domains at each side of u, = 0. We will have then: 


€ d u—e d 
T, (€) = f qo Tr) >) Ta 
x(0) x(u — x?) e x(u — x7) 

(3.13) 

for the y < He and the y > pe; domains, respectively. The 

integrals are performed using different types of initial and final 

x-states. For the first part, where the fixed point is x* =0 we 

start from a given small value x(0) and end up close to x*, 
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Figure 3.9. Divergence of transient times in the symmetry breaking 
problem described by the dynamical system x = ux — x°. Here we 
plot the analytic result obtaind in the text for Te (€) (left curve) and 
Ti) (right), using € = 0.001 and x(0) = 0.05. 


specifically x = € < 1. Similarly, the transient for the second 
half is computed by starting from x(0) = € to yM — €. 

The resulting calculations give our analytic estimates of the 
transient times: 


Lo 1 few = x0) 
ZO= "| Gea) OM 
and 
ol (u — €7)(,/p — €7) 
no-z" |G (Va = I 05 


These two curves are displayed in figure 3.9 against the control 
parameter u. From both sides of the critical point, the transient 
rapidly increases, becoming infinite at ye = 0. For m values very 
close to u, we have in fact T (€) ~ 1/u which gives a hyperbolic 
decay around the critical point. 

The key implication of this behavior is that a dynamical signa- 
ture of criticality is a sharply marked increase in the relaxation 
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time exhibited by systems close to the phase transition. Such 
increase has been reported from different experimental studies, 
from brain (Haken 1996) to ecosystem dynamics (see chapters 
12-13). Actually, it can be used as a warning signal of a close 
critical transition and thus as a potential indicator of future 
changes. 


3.5 Multiple Broken Symmetries 


The bifurcations considered in this book are all related to changes 
leading to different types of behavior. In the real world, they cor- 
respond to transitions between alternative states and emergence 
of new possible attractors, and can occur at very different scales 
(from a cell to a society). But in fact we need to keep in mind 
that in some cases the system undergoes a cascade of bifurcations. 
An example of this situation is biological development (Goodwin 
1994). Through the process of embryo development, millions 
of events happen, from cell divisions to cell death and differen- 
tiation into different cell states. All these events are associated 
with nonlinear dynamical processes where gene-gene and cell- 
cell interactions play a role. This was early envisioned by the 
great biologist Conrad Hal Waddington (Waddington 1957; 
Slack 2003) who thought about the problem of development as 
a multistep series of bifurcations (figure 3.10a). In Waddington’s 
picture, it is possible to visualize the whole process as a ball (the 
embryo’s state) rolling down through a multidimensional space. 
Every fork in this landscape would correspond to a bifurcation. 
What controls the shape of the landscape? Clearly, different 
actors will be involved, including genes and their interactions. 
Such elements would define a complex network both shaping and 
constraining the dynamics of development on the landscape, as 
schematically indicated in figure 3.10b. This picture is a powerful 
one and has been useful in understanding the evolution of novelty 
(Barton et al. 2007). Moreover, it is also an interesting metaphor 
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Figure 3.10. Multiple paths in a complex, nonlinear dynamical 
system, as illustrated in (a) by the so called epigenetic landscape of 
Waddington (see text). Here the underlying dynamical system is the 
process of development of an organism. The state of this process is 
indicated by the position of the ball. The process of development 
would correspond to the ball rolling down through this surface. The 
shape of the landscape would be controlled (b) by the genes involved 
in the process, defining a network of interacting units. 


for other types of complex systems: broken symmetries pervade 
multiple aspects of complexity, including the generation of spatial 
structures, memory, and behavior (Anderson 1972; Palmer 1982; 
Hopfield 1994; Livio 2005). 


4 


PERCOLATION 


4.1 Systems Are Connected 


What makes us label a collection of objects a “system”? The most 
obvious answer is that there is some kind of link between different 
objects such that it makes sense to consider them as part of a larger 
structure. If a given set of species, genes, or individuals contains 
elements that appear somehow isolated from the majority, we will 
probably discard them as part of the global structure. A system is 
thus, roughly speaking, a set of items such that we can find a path 
connecting any two of them. 

Once we agree on the previous minimal requirement, new 
questions emerge. How do the elements forming a given system 
get connected? How does connectivity affect their behavior? 
Answering these questions in a general way provides a broad view 
of the problem. In order to do this, let us consider a simple model 
that illustrates a powerful concept (Sahimi 1994; Stanley et al. 
1999; Stauffer and Aharony 1999; Hinrichsen 2000; Christensen 
and Moloney 2005). Imagine you have an L x L square grid 
where each patch is either empty or occupied by a piece of 
flammable material. This is usually called a “tree” and thus the 
set of trees scattered throughout our lattice can also be dubbed a 
“forest.” Imagine that trees are introduced in any given place with 
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Figure 4.1. Three lattices of 150 x 150 sites are shown, with diff- 
erent densities of occupied sites (here indicated as black squares): (a) 
p =0.55, (b) p = 0.60, (c) p = 0.65. Each site can be occupied 
with the same probability p, and the resulting spatial pattern is 
almost the same except for the effective number of occupied sites. 


some probability p. At low p levels, trees will appear isolated, 
whereas high occupation probabilities will crowd the system, with 
each tree surrounded by several others. In figure 4.1 we show three 
examples of different levels of tree densities. 

Although the density of trees increases linearly with p (we will 
have p L? trees on average) the global connectivity of the system is 
far from linear. Let us define the neighborhood of a given patch as 
the set of four nearest patches. This is the so-called von Neumann 
neighborhood: Trees are connected only to their four neighbors 
and if a tree gets burned, the fire can propagate only to those four 
locations (if trees are present). It seems clear that if a fire starts at 
some given location, it will hardly propagate at low p values and 
instead it will burn most trees at high p’s. Our intuition suggests 
that the probability of a fire spreading over a large area would 
increase somehow linearly with p. Is this the case? 


4.2 Percolation Thresholds 


In order to answer the last question, let us formulate it in more 
quantitative terms. Suppose we burn all trees at the bottom line 
of our lattice and ask what the probability is that the fire will 
reach the upper limit of the lattice. When such an event takes 
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Figure 4.2. Burning trees in the previous figure: here at time zero all 
trees at the bottom line of each plot in the previous figure burn. The 
front propagates and burns further trees. For p = 0.55 (a) the fire 
goes extinct rapidly, whereas for p = 0.65 (c), it spreads to almost 
every tree in the lattice. Close to the critical point (b), a very large 
fire occurs, spanning the whole lattice but leaving many complex 
holes at all scales. The order parameter and the transient times 
associated with the percolation transition are shown in (d) and (e), 
respectively. 


place, we say that the fire percolates the entire system. Percolation 
implies that there is some long path connecting points separated 
by a distance that is of the order of the system’s size. We can 
simulate this process using a computer and calculate the number 
of burned trees for different occupation probabilities. 

The results of these computer experiments are shown in 
figure 4.2. Despite our intuition, we can see that the fraction of 
trees burned is highly nonlinear. The burning actually displays 
a two-phase behavior. Below some threshold (figure 4.2a) the 
fires essentially die out very rapidly. At this subcritical phase 
the low density of trees makes fire propagation too difficult to 
succeed. Instead, once we cross the so-called percolation threshold 
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(figure 4.2b, located at p, ~% 0.59), the fire is very likely to 
percolate the entire lattice, and the size of the fires rapidly grows 
(figure 4.2c). This is rather unexpected and has a number of 
important consequences; the first is that we need to think of 
percolation as a critical phenomenon: it is in fact an example 
of a second-order phase transition. The second is that, in order 
to become a system, the minimal requirements to achieve global 
connectivity involve the presence of a well-defined threshold. 
Once such a threshold is reached, information, fires, or epidemics 
can easily propagate at the large scale. Determining the conditions 
under which percolation occurs is often equivalent to finding the 
boundary between two well-defined phases. Such phases describe 
on one hand decoupled subsystems unable to propagate signals 
and perturbations and, on the other, coherent structures where 
information is exchanged and processed globally. 

This transition can also be described in terms of a control 
and an order parameter, as with the Ising model in chapter 1 
(see figure 1.5a). Here the control parameter is p, whereas the 
order parameter (2(p) can be, for example, the fraction of burned 
sites. Specifically, we run many simulations starting from random 
initial conditions (different trees will burn at the bottom of our 
lattices) and calculate in the end the size of the fire. In figure 4.2d 
we display the result of our numerical simulations for a L = 107 
lattice using 10° runs. The result is familiar to us: a second 
order (critical) transition occurs near p, = 0.59, where fires start 
to successfully propagate. Just increasing p a little beyond this 
point will affect most elements. The sharpness of this transition 
increases dramatically as the system size is increased. For very large 
L, we would see a marked change between an almost-zero-order 
parameter for all p < p, toanon-zero value after the critical point 
is reached. 

Two important comments need to be made. First, the com- 
plexity of the pattern observed at criticality is clearly larger than 
the ones obtained for both p > p, (very homogeneous, random) 


4.3 Percolation Is Widespread 57 


and p < p, (homogeneous, flat). At criticality, multiple scales 
are involved, as illustrated by the presence of islands of patches 
of all sizes, displaying a fractal structure. The complexity of 
this structure can be measured in many ways, but a simple, 
indirect measure is provided by the transient time t(p) required 
in order to complete the burning process. In figure 4.2e the 
largest transient time observed for our simulations is shown. For 
subcritical values of p, the fire dies out very rapidly, as expected, 
whereas for high p values, it is close to L, also consistent with our 
expectations. However, a marked peak appears around the critical 
point: it takes a long time to burn the complex structure shown 
in figure 4.2b. This is due to the complex, nested geometry of the 
percolation cluster. 

The second point is related to the different nature of this 
transition (and most of those considered in this book) in relation 
to the Ising model. The percolation process described above is 
an example of a so-called absorbing phase transition (Hinrichsen 
2000). These types of transitions are characterized by the presence 
of a stationary state that can be reached by the system but from 
which the system cannot escape. This is the most important 
class of nonequilibrium phase transitions, and as we mentioned 
in chapter 1, they violate the detailed balance condition. Such 
a class of behavioral pattern does not allow a connection with 
equilibrium thermodynamics or define energy functions, as in the 
Ising model. 


4.3 Percolation Is Widespread 


The percolation transition appears in many different contexts 
(Stauffer and Aharony 1999; Christensen and Moloney 2005). 
An interesting example is provided by the so-called sol-gel tran- 
sition (Pollack 2001). Let us imagine a given set of complex 
molecules each having f different active centers (also called 
functional groups). These can be seen as different parts of the 
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Figure 4.3. Sol-gel phase transition. In (a) the mean-field parameter 
space is shown, with the two phases (sol/gel) separated by the critical 
line pe =1/(f — 1). In (b) an example of collagen fibers under 
the He-Ion electron microscope is shown (image courtesy of Zeiss 
Nanotechnology). Figure (c) shows a schematic representation of the 
percolation network formed in collagen solutions (Redrawn from 
Forgacs et al. 1991). 


molecule that can react with other molecules, forming a bond 
and thus a larger macromolecular structure. Once a given pair of 
molecules has reacted, each one has f — 1 free centers which can 
also react. If p is the probability that a given bond gets formed, 
what is the expected number of bonds for a given randomly 
chosen molecule? It is easy to see that the average number is 
simply p(f — 1). But now this provides a critical condition for 
percolation to occur: if p(f — 1) < 1 then on average it is likely 
that no further links will be made, whereas if p(f — 1) > 1, then 
at least one link will be expected. The critical boundary is thus 
given by: 
1 
P= Ten 


The resulting separation of the (p, f) parameter space into two 
phases is shown in figure 4.3a. 

This situation appears in a number of different contexts 
involving polymers (Atkins 2000; de Gennes 1976). It has also 
been identified in some key molecular processes associated with 


(4.1) 


morphogenesis (Forgacs et al. 1991; Forgacs and Newman 2007). 


4.4 Bethe Lattices 59 


In some types of tissues (such as mesenchyme)’ the morphological 
changes and structure displayed are largely dependent on the 
properties of the network of fibers that are scattered within 
the extracellular matrix where cells are embedded. A particu- 
larly important component is collagen.” Collagen undergoes an 
assembly process generating macromolecular fibrils and, under 
the appropriate conditions, macroscopic fibers. Forgacs and co- 
workers showed that once a critical density of fibers is reached 
(figure 4.3b), percolation takes place (Forgacs et al. 1991; Forgacs 
1995). 


4.4 Bethe Lattices 


The previous derivation of the percolation point for the sol-gel 
example provides an estimate to the critical conditions under 
which the system becomes connected. What else can be said? A 
significant quantity is the expected size of the largest connected 
set of elements, the so-called giant component (see also chapter 
5). Specifically, although p,(f) gives a necessary condition for 
the existence of a connected system, it is clear that not all elements 
will be connected through some path. But we can calculate the 
fraction of elements belonging to the giant component under 
a mean field approach. This approximation is based on the so- 
called Bethe lattice (Stauffer and Aharony 1992) and provides a 
simple (mean field) model of percolation. This type of lattice has 
a tree structure (figure 4.4a) with no loops that has been used in 
epidemics as a simple illustration of disease propagation. 

Let us consider the nodes of this lattice as occupied or free 
and assume that disease or fire can propagate through from site 


! This is a loosely organized connective tissue displaying a viscous consistency. 
It contains a network of fibers and a class of cells known as fibroblasts. 

? This is the single most abundant molecule in the human body, making up 
as much as 6 percent of our body weight. 


Figure 4.4. Percolation on Bethe lattices. In (a) we show an example 
of the tree structure of a Bethe lattice. Starting from the center 
moving outward, each node has a fixed number of outgoing links. 
In (b) the phase transition, as predicted from mean field theory 
arguments, is shown. 


to site if and only if both sites are both occupied. The branching 
is assumed to be directional, which is to say, although each node 
has z links (the so-called coordination number), only z — 1 of 
them can actually propagate the signal. If p is the probability 
of occupation of a site, the average number of sites that can 
propagate the epidemics will be p(z — 1). The critical condition 
for propagation defining the percolation threshold will be the 
same as before, namely p, (z — 1) = 1. 

Assuming that the lattice is infinite, we can also compute the 
percolation probability P(p) that, given one infected site, the 
disease will propagate (percolate) to infinity. For p < p, we have 
P(p)=0, but for p > pe the specific form of P(p) must be 
derived. If we label as Cœ an infinite (percolating cluster) and 
call Q the probability that a randomly chosen site s; ¢ Ca, it is 
not difficult to see that: 


Q= Pls; ¢ Co] =(1— p) + pQ"' (4.2) 


The first term on the right-hand side corresponds to the proba- 
bility that the node is empty. The second term is the probability 
that s; is occupied and none of the neighbors belongs to Cy. 
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For z = 3, the previous equation gives two solutions: Q = 1 and 


Q= (1 — p)/p. From them the percolation probability will be: 
; lapy 
P(p) = pl — Q=p]|1- Fa (4.3) 


and it is shown in figure 4.4b. A rapid increase in the percolation 
probability is observed at p, = 0.5. Actually, we can find the 
approximated form of P(p) close to p, by using the Taylor 


expansion: 
wey 
(p) p dp), ae 
1 (d’P(p) j 
af Fry ). (=p) pee (44) 


which gives, ignoring higher-order terms, P(p) © 6(p — p-). It 
can be shown that, the larger the value of z, the steeper the 
transition curve. This is consistent with the presence of larger 
numbers of alternative paths able to sustain propagation. 

Other quantities can be derived using apprpriate mean field 
arguments (Christensen and Moloney 2005). As an example, 
consider the problem of finding the average of cluster sizes 
(defined in terms of sets of connected, occupied nodes). It can 
be shown that it is a function x (p) given by 


_ Pc + p) 
Pe 


for p < pe. As we can see, we have again a divergence close to 


x(p) (4.5) 


Pes that is, xX (p) ~ (pe = py! This result is consistent with our 
previous finding: close to percolation thresholds, the likelihood 
of finding a connected path (i.e., chains of connected elements 
within clusters) rapidly grows. 

Percolation and bifurcation analysis can be understood as two 
different forms of detecting and characterizing nonequilibrium 
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phase transitions under mean field approximations. In the next 
chapters, different case studies will be considered and specific an- 
alytic tools will be used, although in many cases both percolation 
and bifurcation theory could be used. 


5 


RANDOM GRAPHS 


5.1 The Erdos-Renyi Model 


As discussed in chapter 4 (on percolation), an interesting phe- 
nomenon that seems to pervade many facets of complexity 
involves the critical requirements that separate a connected from 
a disconnected system. In our examples we used spatially ex- 
tended systems (lattices) showing that, for a threshold level of 
local connectivity, a signal could propagate through the system 
spanning its entire size. A simpler version of this idea does 
not consider space and is thus free from spatial correlations: it 
provides the simplest scenario for percolation within a set of 
connected elements under mean field assumptions. 

Our example is based on a graph Q = (V, E), defined by 
means of a set of N vertices (or nodes) 


V = {vj, v2,..., vy} (5.1) 
and a set of L edges (or links) 
E ={e,,e2,..., ez} (5.2) 


connecting pairs of vertices. For each node v; € V we can define 
its degree k; as the number of edges connecting v; with other 
nodes in the graph. The average degree (k), also indicated as z, 
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is simply 
l 
(ky = k (5.3) 


The simplest model of a random network, known as the Erdos- 
Renyi (ER) graph (Bollobas 1985), allows us to define a null 
model of a disordered network. Although this system is seldom 
found in reality, it provides an excellent case study when looking 
at how networks depart from homogeneity. This type of random 
model has been actually used in different contexts, including eco- 
logical (May 1976; Cohen 1978), genetic or metabolic (Kauffman 
1962, 1993), and neural (Hertz 1991, Sompolinsky 1988, Amari 
1988) networks. 


5.2 Statistical Patterns 


In ER graphs, every edge is present with probability p, and absent 
with probability 1 — p.’ The algorithm for building such web is 
simple: choose every pair of nodes, generate a random number 
E € (0, 1), and add a link if £ < p. 
The average number of edges (L) in these graphs is thus: 
NN — 1) 


He- n f (5.4) 


where the 1/2 factor is due to the fact that each link is shared 
by two edges, and accordingly, the average degree can be written 
as z=2(L)}/ N~ Np being valid only for large N. It can be 
easily shown that the probability P (k) that a vertex has a degree 
k follows a Binomial shape, 


P(k) = (” 5 ') pa pi! (5.5) 


! In fact, there is a whole ensemble of graphs for a single value of p, and random 
graph theory deals with the average properties of these ensembles. 
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Figure 5.1. Random homogeneous networks. In (a) a typical ER 
graph is shown for z= 3. In (b) the corresponding degree distri- 
bution is shown. 


and for large N and fixed z, a good description is provided by a 
Poisson distribution: 
k 


Ph) = e~? 


A (5.6) 


Moreover, the Poisson distribution has an interesting property: 
the variance o(k), defined as of = (k”) — (k), can be shown to 
have the same value as the mean, (k) = øf. These distributions 
are characterized by the average degree z: a randomly chosen 
vertex is likely to have a degree z or close to it. An example of the 
resulting network (for V = 100 and z = 3) is shown in figure 5.1, 
together with the corresponding degree distribution. 


5.3 Percolation at Critical Connectivity 


One of the most interesting properties exhibited by these graphs 
concerns the presence of a percolation threshold, separating a 
disconnected network (where no information can be exchanged 
among components) to a system where most (if not all) elements 
are linked to one another through at least one path. It is clear 
that increasing p in the ER model effectively increases the 


Figure 5.2. Random graphs: four examples of ER graphs are shown, 
for different values of the average degree z = Np with N = 100 (here 
nodes are placed at equal distance along a circle). As p grows, the 
number of isolated elements gets reduced. Such reduction is sharp 
once Z, = 1 is crossed. 


probability of two given nodes of connection through some path. 
Intuitively, for very small p, the graph will be fragmented into 
many components, most pairs being isolated or belonging to a 
pair of connected units. However, if p is large, most pairs of nodes 
will be linked and thus very short paths will be observed. Now the 
question is: as we increase the average number of links z, what is 
the probability of observing a connected (or almost connected) 
graph? Specifically, if S is the fraction of nodes belonging to the 
largest connected component (the so-called giant component), we 
want to know how the function S(p) behaves as p is changed.* 
Against our intuitions, and consistent with our example of a 
forest fire from chapter 4, such probability does not increase in 
a monotonous fashion with z (alternatively, with p, once N is 
fixed, since z = p N). This model displays a phase transition at 
a given critical average degree z, = 1 (figure 5.2). At this critical 
point, it is said that a giant component forms (Bollobas 1985; an 
excellent presentation is given in Newman 2003a): for z > z, a 
large fraction of nodes are connected in a single web, whereas for 


Z < z, the system is fragmented into (many) small subwebs. 


? In all these derivations, we assume the existence of an ensemble Gi, of 
graphs of size N and parameter p. 
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To see how this happens, let us follow a simple approximation 
that shares some points with our previous derivation using the 
Bethe lattice. Consider the graph Q as the joining of a giant 
component S% and the rest of the elements not belonging to 
it. Let us indicate as Q the probability that a randomly chosen 
element v; does not belong to the giant component, that is, 
Q= Pv; ¢ Soo] (figure 5.3). The probability of a given node 
of degree & not belonging to the giant component must be equal 
to the probability that none of its k neighbors belong to Soo, that 
is, Q*. The average value of Q must be estimated by averaging all 
of k, that is, 


Q=X PHQ (5.7) 


[9] 
k=0 


For a Poissonian graph, using the Taylor expansion } `, x/k! = e*, 
we obtain: 


(5.8) 
k=0 

Since the size S of the giant component is S = 1 — Q, we have 

the following transcendental equation for the growth of the giant 

component with z: 


S=1-e* (5.9) 


We can see that for z —> 00 we have S— 1. Besides, we 
can also check that for z= 1 the previous equation is valid for 
S=0. The numerical solution to the last equation is plotted 
in figure 5.3c, showing that the (second order) transition from 
the disconnected (z < 1) to the connected (z > 1) phase occurs 
suddenly at the critical point z,. 

The importance of this phenomenon is obvious in terms of 
the collective properties that arise at the critical point: systems- 
level communication becomes possible and information can flow 
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Figure 5.3. In order to estimate the value of the giant component 
S for the ER graph, we consider the probability of a given node v; 
to belong (a) or not (b) to the giant component. In (c) we show 
the fraction S of the graph belonging to the giant component as 
the average degree z is changed. A phase transition occurs at z; = 1, 
where the fraction of nodes included in the giant component rapidly 
increases. 


from the units to the whole system and back. Moreover, the 
transition occurs suddenly and implies an innovation. No less 
important, it takes place at a low cost in terms of the number 
of required links. Since the only requirement for reaching whole 
communication is to have a given (small) number of links per 
node, once the threshold is reached, order can emerge for free. In 
this context, an additional result from the ER analysis concerns 
global connectedness. Specifically, when p >In N/N (i.e., for 
z >In N), almost any graph in G No is totally connected. 


5.4 Real Networks Are Heterogeneous 


The present example is important in providing a shortcut to 
the percolation problem already discussed in chapter 4. By 
ignoring lattice effects (and thus spatial correlations), we obtain 
an extremely simple model of a system formed by a large number 
of connected units. But we can also ask ourselves how general the 
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ER approach can be. Are real networks at all like this? The answer 
is essentially negative. Real networks have been found to diverge 
from this homogeneous picture in a number of ways (Albert and 
Barabdsi 2002; Dogorovtsev and Mendes 2003; Newman 2003b). 

Perhaps the most interesting feature exhibited by most real 
networks, from the genome to the Internet, is that they are highly 
heterogeneous: they can be described by a degree distribution 
of the form P(k) = (1/Z)k-’ O(R/E) where P(k/E) introduces 
a cut-off at some characteristic scale E£ and Z is a normalization 
constant.? On a log-log display, these distributions will show, for 
high cut-offs, a straight line behavior. As the value of £ grows, fat 
tails develop. 

Three main classes can be defined (Amaral et al. 2000): 
(a) When & is very small, P(k) ~ @(k/é) and thus the link 
distribution is homogeneous. Typically, this would correspond to 
exponential or Gaussian distributions; (b) as € grows, a power law 
with a sharp cut-off is obtained; (c) for large &, scale-free nets are 
observed. The last two cases have been shown to be widespread 
and their topological properties have immediate consequences 
for network robustness and fragility. In particular, since a scale- 
free shape implies that most elements will have just one or a 
few links, the removal of such nodes will have little effect on 
the network (particularly in relation to its connectivity). In this 
context, these webs are expected to be robust against random 
failures. However, the presence of a few, highly connected hubs 
creates some undesirable effects: if they are removed (or fail to 
properly operate) the whole system can experience a collapse. 


3 This normalization is defined as Z = yy ko (£) where K is the 


maximum number of links that a node displays within the network. 


6 


LIFE ORIGINS 


6.1 Prebiotic Evolution 


How life originated on our planet is considered one of the most 
fundamental questions of science (Maynard Smith and Szathmary 
1999; Dyson 1999, Kauffman 1993). The emergence of living 
structures seems inconsistent with the expected patterns of decay 
and increasing disorder displayed by nonliving matter. And yet, 
somehow self-replicating structures appeared and evolved into 
life forms of increasing complexity. Although cells have become 
the basic units of living organization, the initial steps toward 
life likely involved the emergence of some classes of complex 
chemical reactions leading to selection for efficient replicators not 
necessarily associated with closed membranes (Maynard Smith 
and Szathmary 1999; Schuster 1996; Ricardo and Szostak 2009). 

Under prebiotic conditions, nonliving matter was composed 
of a rich soup of diverse molecules, some of them generated 
on the earth’s surface and others—as suggested by the Catalan 
astrochemist Joan Ord—likely coming from comets falling on 
our planet (Oró 1961a). Such molecules included the building 
blocks of life, for instance, amino acids, nucleotides, lipids, and 
sugars, most of which have been shown to spontaneously emerge 
under laboratory conditions out of water, methane, hydrogen, 
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and other nonbiological molecules (Miller 1953, Oró 1961b). 
But such precursors are not enough for life to emerge: in order to 
go beyond the soup, some chemical reactions need to be amplified 
so that the abundance of certain special molecules increases and 
allows selection processes to be at play (Eigen 1999). If such 
amplification occurs, growth of a subset of chemical species would 
take place. And growth is clearly a first ingredient necessary for 
reproduction. Afterward, Darwinian selection can come into play 
(Szathmáry 2006). 

Growth is easily achieved through autocatalysis when a given 
component A can make a copy of itself by using available 
monomers Æ, that is, 


A+E524+W (6.1) 


being W waste material and u the reaction rate. This system 
can be described using a Malthusian growth process, that is, if 
x = [A], we have dx/dt= ux which gives (see chapter 2) an 
exponential growth dynamics, which is to say, x(t) = xoe”* which 
is obviously valid only for some restricted conditions, since sooner 
or later spatial constraints will slow the dynamics (such as in a 
logistic system). 


6.2 Replicators and Reproducers 


Once autocatalysis is at work, the question is what type of 
complex chemical structures can be sustained and what are the 
minimal requirements for such autocatalytic systems to emerge. 
Cooperation has been suggested as an essential step toward the 
emergence of complex, self-organized chemical systems (Eigen 
1971). In particular, we can consider the effects of having 
cooperation among molecules in order to generate new copies, 
namely reactions such as 


SALES 3A4A+W (6.2) 
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For this system, it is not difficult to see that the corresponding 
equation is dx /dt = ux? and the solution now takes the form 
1 


= 1/xo — ut 


It can be easily shown that this solution has a well-defined 


x(t) (6.3) 


asymptote at a finite time 4, = 1/uxo where a divergence in 
the concentration is observed. This means that after a well- 
defined time interval, an infinite concentration will be reached. 
This type of dynamics is known as /yperbolic and as we can 
see is dependent upon the initial state. Through competi- 
tion among several chemical species (or replicators), the ones 
with higher initial populations would be more likely to win 
the race. 

Here we want to analyze the properties displayed by a system 
where replicators (either molecules or cells) show both Malthusian 
and cooperative dynamics, following the work by Ferreira and 
Fontanari (2002). 


6.3 Autocatalysis and Cooperation 


We will consider a simple situation (figure 6.1) where the repli- 
cators A change their concentrations through autocatalysis, coop- 
eration, and decay: 


A= 2A 2434 A350 (6.4) 


where the last one occurs at a unit rate. For simplicity, we do 
not include the precursors or waste materials here (although both 
are always present) since they would not affect our calculations 
provided that precursors are abundant. 

In the following, we will assume that reactions occur in a 
system with limited space and that the available monomers are 
present in high numbers. For convenience we also fix [F] = 1. 
If we consider a normalized concentration 0 < x < 1 using the 


Figure 6.1. A simple model of replicating molecules considers three 
basic reactions, here indicating as happening on a two-dimensional 
lattice. Reactions include: (a) decay, where a given particle (gray 
sphere) spontaneously degrades (empty circle); (b) autocatalysis, 
where a given particle occupies an empty spot and (c) cooperative 
replication, where the presence of a molecule in the neighborhood 
helps the replication of another molecule. 


previous set of reactions, we have a mean field model for x = [A] 
dx 


— = —x +5x(1 — x) +px?(1 — x) (6.5) 
dt 


where the three terms in the right-hand side of the equation cor- 
respond to the three mechanisms indicated in the previous figure. 
An obvious stationary state corresponds to the extinction con- 
figuration (i.e., x* = 0). The stability of this state is obtained from 
"e = =—145(1—2x)+px(2—3x) (6.6) 
which gives 4,,(0) =s — 1, and thus the extinction state will be 
stable ifs < 1. 
Two additional fixed points can be present, and correspond 
to the solutions of the quadratic equation —1 + s(1 — x) + 
pux(1 — x) = 0. This gives: 


1 
saz e-s t Vut] (6.7) 
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By analyzing this expression, it is not difficult to find the 
conditions for the presence or absence of the nontrivial fixed 
points. First, real solutions will exist provided that 


(u +s)? > 4 (6.8) 


and these solutions will be positive (and thus meaningful) pro- 
vided that the inequality u — s > 0 holds, and thus an additional 
constraint to consider is y > s. From the condition of real roots, 
we can find the critical boundary: 


6) = =p 2h (6.9) 


These conditions allow defining three phases in the (u, s) 
parameter space, indicated in figure 6.2a. These correspond to: 
(a) an active phase where molecules are always present (upper 
domain, light gray), (b) an empty phase where the population 
goes extinct, and (c) a bistable phase (dark gray) where two 
alternate states are possible. The last two phases are separated 
by the critical curve s (u) defining the boundary for a first-order 
transition. 

The bistable character of our model is clearly observed in the 
bifurcation diagram shown in figure 6.2b, where we plot the 
fixed points against the cooperative replication rate u for a fixed 
s value (here s = 0.75). As predicted by the stability analysis, 
the extinction point always appears stable: small population 
deviations will fall off to zero. A marked change occurs at a critical 
rate, which can be obtained from our previous equations for each 
s value and is simply u, = 2 — s + 2/71 =s. 

The presence of a first-order transition has two important 
implications for our problem. First, either an efficient autocataly- 
sis or a combination of autocatalytic and cooperative effects needs 
to be present in order to sustain a stable population. Second, a 
minimal initial population must be stabilized in order to avoid 
falling into the extinct state. We can apreciate this situation by 
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Figure 6.2. Phases in the model of molecular replicators. In (a) the 
three basic phases are shown. The black line separating the empty 
(extinction) phase from the bistable phase is given by equation 6.9. 
In (b) an example of the bifurcation diagram is shown for s = 0.75. 
For this parameter value, we have a bifurcation at u, = 0.25. 


using the potential function associated with our system, namely: 


2 3 4 


O,) =U) > +6 -wW> +T (6.10) 


which is plotted in figure 6.3. Here we use s =0.75 and three 
different values of u. The minima defining the alive phase 
coexist with an alternative minimum where extinction is also an 
alternative possibility. When u < u, = 2.25, a unique minimum 
is observable, associated with the extinction scenario (or dead 
phase), whereas for u > ue, we will observe two minima, being 
the alive fixed point placed in a deeper valley. Moreover, since 
fluctuations are expected to be involved in any of these prebiotic 
scenarios, noise could push the system from the alive to the dead 
state, thus making the transition less robust than predicted from 
our deterministic framework. 


6.4 Prebiotic Reaction Networks 


The previous results, as well as others suggested by theoretical 
models (see next chapter) indicate that early life might have 
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Figure 6.3. Stable states in the replicator model describe two basic 
types of prebiotic organization, namely alive and dead. The first 
corresponds to a stable, nonzero concentration of A molecules, 
whereas the second deals with extinction. If the reactivity is too low, 
decay will win over synthesis, but even when the reaction rates are 
appropriate to sustain a stable population, a minimal initial number 
is required. The potential function ®,,(x) provides a good picture of 
the problem. Here three examples are shown (see text). 


faced some serious obstacles. The need for a minimal amount 
of molecules allowing the system to reach a high concentration 
of self-sustaining molecules might indicate that the jump from 
nonliving to living matter was difficult. However, it is known 
that life originated very early on our planet, almost as soon as 
its surface cooled enough to allow a stable chemistry to exist. 

In trying to understand the possible origins of life, some 
theories have used the idea that molecular species can interact in 
diverse ways and create cooperative structures (Eigen and Schuster 
1979; Farmer et al. 1986; Kauffman 1993). The underlying idea 
is that a molecular soup receiving external energy inputs and 
displaying enough chemical diversity will eventually generate a 
self-sustained autocatalytic set. This argument is based on two 
key assumptions. First, available molecules must exhibit catalytic 
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properties, that is, they must actively contribute to allow for new 
reactions to occur. The second is that, as the molecular size of 
the polymers grows, the number of chains that can be created 
explodes. In his early contribution to this problem, Kauffman 
derived a percolation condition for these autocatalytic sets to 
emerge (Kauffman 1993). Once such a percolation threshold 
was reached, complex networks would spontaneously form and 
persist. These networks would also sustain a high diversity of 
chemical molecules. Such diversity can be properly analyzed by 
using the phase transition approach (Jain and Krishna 2001; 
Hanel et al. 2005; see also Szathmáry 2006). 


T 


VIRUS DYNAMICS 


7.1 Molecular Parasites 


Life is rich in variety. Over hundreds of millions of years, life 
forms have been emerging in all habitats around the world. At all 
scales, from nanometers to meters, the success of life is obvious 
and overwhelming. Complexity has been achieved by means 
of different sources of change and innovation. These include 
multicellularity, cooperation, and language. But there is a dark 
side common to all these life forms: the presence of parasites. 
Parasitic entities seem a universal outcome of life. No known 
living species seems to be free from at least one type of parasite. 
The most commons forms of these entities are viruses, which are 
molecular machines that exploit cellular complexity to replicate 
themselves. Their genomes are typically small, too small in any 
case to allow them autonomy, but large enough to allow them to 
enter their hosts and use the host’s molecular machinery for their 
own purposes. 

Viruses have a large diversity of structures and genomes 
(figure 7.1). In terms of genetic material, there are DNA and 
RNA viruses with single-stranded, double-stranded, and frag- 
mented genomes. RNA viruses offer a unique opportunity for 
exploring long-term evolution under controlled conditions due 
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Figure 7.1. In (a) we show an example of RNA virus (HIV) 
where different proteins and other components are indicated. RNA 
viruses display high mutation rates, which scale as the inverse 
of their genome length. This pattern of decay is roughly fol- 
lowed by other organisms, as shown in (b) where we plot the 
logarithm of mutation rate u against the logarithm of genome 
size. The dotted line is the predicted inverse relation among 
both quantities (see text). Figure (a) has been adapted from 
commons.wikimedia.org/wiki/VIHsanslibel.png. 


to their high mutation rates. They are the most important class 
of intracellular parasites and are responsible for a wide range 
of diseases (Domingo 2006). They are extremely adaptable to 
changing conditions, particularly to selection barriers imposed by 
immune responses. The immune system is able to identify pieces 
of the infecting virus, eventually wiping out the original strain. 
However, if the virus changes fast enough, it can escape from 
immune recognition. Since high mutation rates allow viruses to 
escape, one might conclude that the higher the mutation, the 
better. Actually, DNA-based life forms are known to support 
sophisticated mechanisms for error correction, whereas RNA 
viruses lack such molecular machinery. 

The enzymes allowing RNA viral genomes to replicate are in 
fact rather inaccurate. What are the limits to such error rates? We 
know that mutations have deleterious effects: a single change in a 
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nucleotide can easily make the new sequence unable to replicate 
at all. Does such a conflict between change and conservation 
exhibit some optimal trade-off? In the early 1970s, a theoretical 
model of replicating and mutating populations revealed that there 
must be an error threshold beyond which information cannot 
be preserved. It was suggested that due to mutation, these pop- 
ulations or guasispecies would be extremely heterogenous (Eigen 
1971; Mas et al. 2010). The quasispecies structure has numerous 
implications for the biology and associated pathology of RNA 
viruses, the most relevant of which is that the heterogeneous 
population structure is a reservoir of variants with potentially 
useful phenotypes in the face of environmental change. 

Mutation rates cannot reach arbitrary values: there is a thresh- 
old to mutational change beyond which no selection is possible. 
This threshold has been dubbed the error threshold (Schuster 
1994; Domingo 2006). Roughly, it was predicted that the mu- 
tation rate Uy, per nucleotide and replication round should be 
us ~ 1/v, with v the sequence length. The Eigen-Schuster 
theory predicts that, beyond u., genomic information is lost as 
the population enters into a drift phase (Eigen 1971; Schuster 
1994). Such prediction received strong support from the first 
estimations of viral mutation rates for different RNA viruses. 
Actually, it turns out that RNA viruses live right at the error 
threshold. In this chapter we will explain the origins of this 
behavior within the context of phase transitions. As will be 
shown, the boundary separating information preservation and 
randomness is a sharp one. 


7.2 Exploring the Hypercube 


Let us first consider a simulation experiment that provides some 
intuitions about the error catastrophe and its origins. Modeling 
virus populations would require considering cells, how they 
recognize their host, how they enter the cell, how their molecules 
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are copied, how they exchange information with the host cell, and 
how the new, mutated copies, leave the cell. However, we might 
not want to know all the details of their life cycle. Instead, we 
ascertain the limits of the mutability of a given RNA virus and 
how these relate to genome size. As usual, we need to define a 
toy model, where a simple representation of the virus genome is 
taken into account and a complete abstraction of the underlying 
host properties is considered. In such description, a population of 
viruses is considered. Viruses will be strings of v units and will 
replicate with a given accuracy. 

In our approximation, a digitial virus is defined as a string S; 
with k = 1,..., N, being N the total population size. Each string 
represents a small genome of size v, that is, 


Sa), 2.45% (7.1) 


where each element is a Boolean variable, Si, € {0, 1}. Formally, 
we can think of each possible string as a vertex S€ H” of a 
v-dimensional hypercube (see figure 7.2). The hypercube is an 
extremely large object. The examples shown in figure 7.2 illustrate 
how we can construct such a landscape under our simplistic 
assumptions. But for a real virus, even the smallest ones, v is a 
large number, ranging from v * 400 for plant viroids (see Gago 
et al. 2009 and references therein) to v ~ 104 for most well-known 
RNA viruses responsible for infectious diseases.’ 

Let us allow the system to change and evolve. Our model starts 
with an initial population of master sequences and repeats, at each 
generation, JV times, the following set of rules: 


1. We take a string S; at random from the population and 
replicate it with probability f(S;). Here two replication 
probabilities are also defined, one for the master 


' This means that the number of nodes in H” will be hyperastronomic, much 
larger in fact than the number of atoms in our universe. 
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Figure 7.2. Quasispecies live in a hypercube. In its simplest form, 
strings are represented as chains of bits of fixed length v. Here the 
hypercubes H” for v = 1, 2, 3, 4 bits are shown. Populations flow 
between the vertices of each hypercube, being two vertices connected 
through single point mutations. 


sequence (where S% = l forall $ = 1,..., v) and the 
other for the rest of the strings, to be indicated as f, and 
fo respectively. Hereafter we indicate the master 
sequence as € = (1, 1,..., 1) and assume that fi > f. 

2. Replication takes place by replacing one of the strings in 
the population (also chosen at random), say S; 4 S; by 
a copy of S;. The copy mechanism takes place with 
some errors. Such errors take place with some 
probability (or mutation rate 14) per bit and replication 
cycle respectively. In this way, each bit S* can be 
replaced by 1 — S% with probability ju. 


The result of our simulations is shown in figure 7.3, where 
we plot three examples of the time evolution of the number of 
strings, starting from our special initial condition. A different 
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Figure 7.3. Time evolution of the bit string quasispecies model 
described in the text. Here we use a population of N = 500 strings, 
with fi = l and f = 0.25. For small mutation rates (left plot) the 
master sequence (thick line) dominates the population, followed by 
the populations of strings one-bit distant in sequence space. These 
are lumped together in a single population, here indicated as a thin 
continuous line. Other populations with strings differing two, three, 
or more bits from the master are also indicated using dotted lines. As 
we increase u, the master sequence becomes less and less frequent: 
in the middle plot we can see that it falls below other populations, 
but still persists. For higher mutation rates, the master sequence 
disappears (right panel). 


representation is provided by the stationary values of each 
population (figure 7.4). For low mutation rates, the master se- 
quence remains high (thick line) whereas other mutant sequences 
have lower populations. For simplicity, we aggregate all strings 
differing one, two, or more bits from the master string. This is 
typically measured by means of the Hamming distance dy: 


du(S;,S;) = Y |S? - Sil (7.2) 
k=1 


As u, grows, we can observe that the master sequence starts to fall 
below that of other populations, eventually disappearing beyond 
some critical value. 
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Figure 7.4. The error catastrophe: using a bit string model, with 
N = 100 strings of size v = 16, starting from an initial condi- 
tion involving only master sequences. After 2000 x WN replications 
with mutation, we plot the number of master sequences (thick 
line) as well as the number of strings with a Hamming distance 
dy =1,2,..., 16 from the master. Although the master dominates 
at low mutation rates, a sharp change takes place at a critical 
mutation rate, where the master sequence vanishes as the quasispecies 
starts wandering through sequence space at random. Here linear (a) 
and linear-logarithmic (b) plots have been used. 


7.3 Swetina-Schuster Quasispecies Model 


In its simplest form, we can consider a reduced system of 
equations defining a population formed by two basic groups: the 
master sequence xı and the other sequences, which we assume 
to be grouped into an “average” sequence with population x2 
(Schuster 1994). Let us also assume (as a first approximation) 
that mutations occur from the master to the second compartment 
but not in the reverse sense. The enormous size of the sequence 
space H” makes this assumption a good first approximation. 
Each time a master sequence replicates, there are v possible 
neighboring mutants. If we consider a real virus, this means 
a very large number of ways of leaving the master node and 
diffusing through sequence space. To move out from a given node 
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Figure 7.5. Phase transitions in the Swetina-Schuster model. In (a) 
we summarize the basic structure of the two-compartment model. 
In (b) the two phases are shown, as obtained from the mean field 
approximation. Here we represent mutation rate per genome (i.e., 
the probability that some unit in the genome gets mutated) against 
replication rate of the master sequence. The gray region indicates 
the disordered (drift) phase, where the virus population wanders 
randomly through sequence space. The white domain indicates the 
presence of a stable cloud of strings around the master sequence. 


is thus easy, whereas going back is unlikely. The basic scheme 
is summarized in figure 7.5a. The two populations are linked 
through mutation from xı to xz (but not back). Additionally, 
each population replicates at well-defined rates, and an out- 
flow rate ® is also included in order to maintain populations at 
a finite size. 

The model is given by the following set of coupled differential 


equations: 


dx, 

e fil — w)x1 — x1®(x1, x2) (7.3) 
d. 

<= fiumi + fiamb) (74) 


where y is the mutation rate, fı is the master sequence replication 
rate, and f is the replication rate for the other strings. The 


86 Virus Dynamics 


term (x1, x2) will allow selection to be effective. This is done 
by assuming a constant population (CP) constraint, namely that 
xı + x2 is kept constant. This means that 


d(x + x2) _ dx dx 


=0 : 
dt dt dt ee 
which allows defining the outflow term as: 
®(n,x,) = Lat he 7.6) 
xy + x2 


Since the relative concentrations of x; and x are given by 


x1 x2 
P = 


P = 2 
xı + x2 xi + x2 


(7.7) 


we have in fact demonstrated that the outflow is nothing but the 
average replication rate: 


(x, x) = Pf tehh = (f) (7.8) 


For simplicity, we fix xı + x2 = 1 in such a way that the ~x;’s are 
already probabilities. Using this constraint, the previous model 
can be reduced to a one-dimensional system (using x2 = 1 — x1) 
and this leads to: 
dx _ 
dt 
where €) =1 — u — fo/fi and &=1 — fo/fi. Here the fixed 


points are xò = 0, xf = &/&, and thus the nontrivial fixed point, 


fixilé: — &x1] (7.9) 


representing a nonzero master sequence population, will be: 


x=] Lf 
fi-f 


The linear stability analysis shows that x = 0 is stable (and 


(7.10) 


the master sequence disappears) if 


b> be =1l- = (7.11) 
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whereas a stable cloud of mutants coexisting with the master 
sequence will be stable for u < ue. The critical mutation rate 
separates the two basic phases. In figure 7.5b we show these 
phases, separated by the critical line, using a fixed fọ value 
(here we use f) = 0.25). Once such a boundary is crossed, we 
shift from one type of qualitative dynamics to the other. The 
standard-error threshold condition is associated with an increase 
in mutation rate. Increased u values crossing the critical line drive 
the master sequence into extinction. 


7.4 Critical Genome Size 


A different form of the previous critical condition can be derived 
by explicitly using the v parameter. This derivation will allow 
us to recover the inverse relation between mutation rates and 
genome size. In order to obtain such a result, let us first write jx as 
a function of v and u, (mutation rate per bit). It is not difficult 
to see that they are related through: 


w=1—(1— u)” (7.12) 


where we have used the probability of no mutation of a given 
unit, 1 — u+, and considered the probability p = (1 — u)” that 
none of the nucleotides is mutated. The difference 1 — p is just 
the probability that some unit does mutate. Since u, is small, we 
can use the approximation’: 


wrl—e MY (7.13) 


It is not difficult to show that the previous critical condition is 
now written, for the mutation rate per unit, as 
a 


hy, = — (7.14) 
Y 


? Here we make use of the first two terms of the Taylor expansion of the 
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Figure 7.6. The two phases displayed by the Swetina-Schuster 
model using the mutation rate per unit and round of replication 
against gemome size. The critical line separating the two phases in 
(a) gives the mutation-genome size observed in real RNA viruses, 
shown in (b) in log-log scale. The resulting power law shows that the 
critical mutation scales as the inverse of genome length. 


where œ = — In(f2/fi)>0 is a constant. The last expression 
actually corresponds to the observed inverse decay of mutation 
rates in RNA viruses as an inverse of their genome size. The 
resulting parameter space is also shown in figure 7.6a. Moreover, 
if we plot the critical line u4(v) in log-log scale (here we use 
fi = 1.0), a well-defined power law is recovered (figure 7.6b). An 
alternate derivation of this condition using a simple percolation- 
like argument can be obtained based on branching processes 
(Demetrius et al. 1985; Hofbauer and Sigmund 1991). 


7.5 Beyond the Threshold 


Increased mutagenesis beyond the error catastrophe can destroy 
the virus, since beyond the threshold no Darwinian selection is 
at work (Schuster 1994). Moreover, considering the fact that 
many mutants will be lethal (at any mutation rate), several critical 


thresholds can be obtained combining lethality and catastrophe 
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(Bull et al. 2005). The exceptionally high mutation rate in RNA 
viruses is illustrated by the finding that most HIV virions in 
blood appear to be nonviable. Effective experimental strategies 
have shown that the error threshold can actually be exploited 
in antiviral therapy (Loeb et al. 2003). Within the context of 
HIV treatment, using promutagenic nucleoside analogs, viral 
replication of HIV has been shown to be abolished in vitro 
(Loeb et al. 2003). Moreover, it has been suggested that genetic 
instability in cancer cells might lead to an error catastrophe: 
mutations and chromosomal changes can increase up to some 
limit (Solé 2003; Solé and Deisboeck 2004; Gatenby and Frieden 
2002; see also chapter 11). 


Our previous models consider a simplistic fitness landscape where 
one given string has a large fitness in relation to all other possible 
strings. Such situation is of course an oversimplification. Fitness 
landscapes are in general much more complex and rich (Bull 
et al. 2005; Jain and Krug 2007), and other components such 
as space (Altmeyer and McCaskill 2001; Pastor-Satorras and Solé 
2001) and recombination (Boerjlist et al. 1996) can also play an 
important role. 

Many developments of the quasispecies model have been 
considered since its original formulation, considering the effects 
of multiple peaks, correlations, and time-dependent changes. In 
this context, an important finding was the presence of neutrality 
in RNA landscapes (van Nimwegen et al. 1999; Wilke 2001a, 
b). In a nutshell, many mutations are selectively neutral (Kimura 
1983) and it has also been found that chains of single muta- 
tions allow connections between neutral genotypes in sequence 
space (our hypercube) in such a way that these neutral net- 
works percolate through large domains of such space (Schuster 
et al. 1994; van Nimweggen et al. 1999; van Nimweggen and 
Crutchfield 2000; Reydis et al. 2001). Neutral networks allow 
populations to display high robustness against mutations and have 
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some unexpected implications. One of them is that, under high 
mutation rates, selection can sometimes favor slower—but more 
robust—genomes. Such a survival of the flattest effect has been 
found both in simulated (Wilke et al. 2001; Sardanyes et al. 2008) 
and real (Codofer et al. 2006; Sanjuan et al. 2007) systems. 


8 


CELL STRUCTURE 


8.1 The Architecture of Cells 


Cells are the basic units of life. They have often been compared 
to machines, involving many components in interaction. These 
components, particularly proteins, define a set of functional 
structures able to cope with all cell requirements, including the 
gathering and processing of matter, energy, and information 
(Harold 2001). Moreover, they also sustain the internal structural 
organization of cells: proteins are the basic units from which com- 
plex chains of polymers are formed. Some of these polymers allow 
cells to define their shapes. They form the so-called cytoskeleton 
(Alberts et al. 2004; Pollard 2003). This skeleton functions a 
scaffold, a positioning system, a distribution grid, and a force- 
generating apparatus used for locomotion (Karp 1999; Pollack 
2001). Most of this molecular machinery is constantly changing, 
and its behavior lies somewhere between order and disorder. 
Understanding how they coordinately organize in space and time 
is a great challenge to experimental and theoretical approaches 


alike. 
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Figure 8.1. Microtubules are formed by tubulin dimers (a) involv- 
ing two almost identical subunits, labeled œ and £ respectively. 
These dimers polymerize, eventually forming long, cylindrical struc- 
tures. In (b-c) the dynamical patterns associated with growth and 
decay of the microtubule are shown (adapted from Howard and 
Hyman 2003). They grow by addition of GTP-bound tubuline, 
and shrink through hydrolysis. If hydrolysis dominates, the cap 
disappears and the structure rapidly dissociates. 


Microtubules are one particularly important component of cell 
organization.' They are not essential only to cell architecture. 
They are highly dynamical, in that they are involved in several 
key cellular processes such as cell division, motion, and transport 
(Howard and Hyman 2003). They are an essential factor for 
cell stability, and actually some anticancer treatments involve 
molecular interactions between drugs and the proteins forming 
the microtubules. They can be properly described as cell railways, 
since they provide the physical substrate along which transport 
of cargo takes place. The basic molecular organization of these 
complex polymers is outlined in figure 8.1. The dynamical 


'The other two are the actin filaments and the so-called intermediate 
filaments. The two involve a diverse range of biological polymers. All three 
groups interact through additional proteins. 
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properties of microtubule formation strongly differ from other 
self-assembling systems, such as ribosomes and viruses. The 
latter involve stable structures, whereas microtubules are in a 
constant state of change and display a broad range of lengths. As 
pointed out in (Gerhardt and Kirshner 1997) they seem to avoid 
equilibrium. 

Microtubules are long, hollow tubular structures. They are 
assembled from dimeric building blocks, combining two very 
similar proteins (figure 8.1a), so-called @-tubulin and 6-tubulin, 
respectively. Each microtubule is like a cylinder composed of 
thirteen parallel filaments having a well-defined polarity. The 
two ends of the tube are usually indicated as the + and — 
ends. The elongation of the microtubules takes place at the + 
end. Polymerization and depolymerization take place all the time 
within cells, and this process displays many fluctuations. A given 
microtubule can grow steadily over a long period of time and then 
rapidly shrink, a behavior known as dynamic instability. 

The driving force of the dynamic instability is the transforma- 
tion of a small molecule, G TP, into another molecule GDP.” 
GTP binds to both tubulin subunits, but only those associated 
to the $ can experience the GTP — GDP change, also known 
as the hydrolyzation of G TP. When the GTP is hydrolyzed, the 
chain becomes curved and thus the cylinder becomes less stable. 
This process is delayed and can take place at different speeds 
depending on local conditions. If the polymerization process is 
faster than the rate of GTP — GDP conversion, growth will 
occur. If not, the microtubule will shrink. The balance between 
these two basic processes defines two possible phases. A cell can 
tune these rates accordingly in order to favor growth or allow 
shape changes. Our goal here is to find the boundary separating 
these two phases using a minimal model of polymerization 
dynamics. 


? GTP and GDP stand for guanidine tri- and di-phosphate, respectively. 
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8.2 States and Transitions 


In order to define a model of microtubule dynamics that could be 
solved in terms of a simple, one-dimensional dynamical system, 
several drastic assumptions must be made. Here we follow the 
approach taken by Antal, Krapivsky, and Redner (AKR) to 
tackle this problem (Antal et al. 2007a, 2007b). First, we will 
consider our problem as an essentially one-dimensional one, 
where a microtubule would be a linear string. Moreover, instead 
of worrying about the dimeric nature of tubulin, we will simply 
consider each bead in the linear chain as a monomer. This 
decision seems reasonable since the relevant component affects 
just one of the tubulins (the £ unit). 

Using this view, we now concentrate on what is going on at the 
tip of the string. In particular, we can look at the end of the string 
and see what is located in the last position. We will indicate as 
filled and empty circles the GTP and GDP subunits, respectively. 
Following the AKR approach, we will indicate the two possible 
terminal states as follows: 


|---@) (8.1) 
for the GTP ending tip and 
jeso) (8.2) 


for the GDP one. The dots simply indicate the existence of a 
linear chain, no matter its composition. Additionally, we will use 
|...) to indicate any string, independent of how its tip is defined. 

Several possible key transitions can be defined here, associated 
to how the tip of our microtubule changes. In the AKR model, 
the following four changes were considered. 


1. A GTP end can incorporate a new GTP at a rate A: 


(or er (8.3) 
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Figure 8.2. Microtubules are constantly changing. Two basic types 
of dynamical behaviors are at work: either they grow smoothly by 
adding new tubulin dimers (a) or rapidly shrink (b). 


2. A GDP end can grow by incorporating a GTP and at a 
rate pA: 


ees ae eee (8.4) 


3. Any intermediate (nonending) GTP unit can transform 
into a GDP end at a rate one: 


E E E [recead (8.5) 
4. A GDP ending can be removed at a rate p: 
ES E (8.6) 
This set of rules can be used as the basis for a mean field 


approximation, where the details of microtubule structure and 
fluctuations are ignored. 


8.3 Dynamical Instability Model 


Two key quantities need to be taken into account in order to 
derive the mean field model for this process. First, the fraction of 
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GDP ends, that is, the probability of having a o at the tip, which 


we will indicate as 
no = P||---o)] (8.7) 


By using the previous list of transitions, we can write a rate 
equation for the dynamics of no. This equation reads: 


d 
ZR = —pàm + (1 — mo) — UN (8.8) 
t 
where we indicate as 
No =|-:-@0) (8.9) 


that is, the frequency of terminal tips ending in a GTP-GDP pair 
(in this order). 

The second component in our analytic derivation is the speed 
V(A, u, p) at which the tip of the string is growing: 


V(A, u, p) = pang + AC — no) — uno (8.10) 


The differential equation for no cannot be easily solved, since 
the last term Mo would be difficult to determine. However, we can 
make a safe approximation by considering that the confirguration 
| +++ @ o) is uncommon, since the GTP units not located at the 
end of the microtubule are rapidly transferred and thus No ~ 0. 
With this approximation, our previous dynamical equation now 


reads: 
dn 
aE = 1 — (1 + pa)no (8.11) 
which gives a single and stable fixed point: 
1 
i= (8.12) 
I+ pa 


Now our interest is to find the critical conditions separating 
the two dynamical phases. Here the speed of growth equation 
will be very useful in defining a percolation-like condition: the tip 
will grow (or shrink) provided that V(A, u, p) is larger (smaller) 
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Figure 8.3. Phases of microtubule behavior from the mean field 
model. The two main types of behavior are described by the two areas 
separated by the critical curve u, = pA(1 + A) (here for p = 1). 


than zero. The critical set of parameters is thus obtained from: 
pano + AC — no) — uno = 0 (8.13) 


and now we can use the stationary state 75, as previously calcu- 
lated. The resulting critical curve is: 


Ite = pM +A) (8.14) 


When u >u. we have the compact phase, whereas the grow 
phase is observed for u< m.. These phases are shown in 
figure 8.3. Although the assumptions used above are rather 
strong, it can be shown that these results are a good approxi- 
mation of the exact behavior (Antal et al. 2007a, 2007b). 


8.4 Tensegrity 


As mentioned at the beginning of this chapter, microtubules 
are one component of a large family of polymers involved in 
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shaping cellular architecture. Most of them are highly dynamical, 
although different time scales are actually involved. They define 
a limited range of molecular mechanisms controlling cellular 
changes (Pollack 2001). The observation that a cell’s cytoplasm 
behaves very much as a polymer makes it possible to explore its 
adaptive behavior and plasticity in terms of simple models. 

Beyond microtubule dynamics, several relevant problems 
remain open within cell biology, particularly in relation to large- 
scale patterns of cell organization. One approach to this problem 
has been made by means of so-called tensegrity, a building 
principle originally described with in the context of architecture 
(Fuller 1961). Tensegrity systems can be defined as structures that 
stabilize their shape through continuous tension. The concept 
has been applied to cell biology (Ingber 1998; 2003) and allows 
connecting mechanical forces and biochemical processes, thus 
integrating global behavior and cell regulation. Given the key role 
played by dynamical transitions between different phases in this 
type of polymer (such as sol-gel changes, see Pollack 2001), phase 
transition models should play a key role in our understanding of 
cell architecture and its dynamical changes. 


9 


EPIDEMIC SPREADING 


9.1 Spreading Diseases 


Infectious diseases have played an important role throughout the 
entire history of human life (Garrett 1994; Diamond 1999). They 
have been a significant source of mortality and have served as 
strong selective forces. No example better summarizes the terrible 
effect of infectious diseases than the Black Plague, which caused 
terror throughout Europe during the years 1347—50. Upon its 
arrival from the East, it spread through Europe killing about a 
quarter of the population. The disease was terribly virulent and 
effective. About 80 percent of the people who were in contact 
with the disease died within two to three days. 

Infectious diseases endemic in Europe such as measles, small- 
pox, influenza, and bubonic plague were transmitted by invad- 
ing peoples and have decimated entire ethnic groups that had 
not been in contact with the diseases and so had evolved no 
immunity. These infectious diseases played an influential role in 
European conquests (Diamond 1999). For example, a smallpox 
epidemic devastated the Aztecs after the first attack by the 
Spaniards in 1520. It is believed that up to 95 percent of the pre- 
Columbian Native American population was killed by diseases 
introduced by Europeans. And similar outcomes took place in 
other parts of the globe. 
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Epidemic modeling has been a very active area of research for 
several decades. One of the simplest model approaches to this 
problem is the so-called contact process. In this simple frame- 
work, two types of populations are considered, namely infective 
(I) and susceptible (S) individuals. The first are infected and carry 
the pathogen with them while the second are not infected but can 
become so. The states of these individuals can change due to two 
basic processes: recovery (J —> S) and infection (S — J), which 
occur with some fixedrates.' In figure 9.1a, the basic transitions 
are summarized for a one-dimensional case (a realistic situation 
requires at least two dimensions). While the recovery process is 
independent of the state of the neighbors, the infection process is 
not: the more neighbors are infected, the larger the probability of 
infection. 

In figure 9.1b we show the results of a numerical simulation 
where we start from an initial condition with half (randomly 
chosen) individuals being infected (here indicated as black dots). 
Here we use V = 100 individuals. Three examples of the dy- 
namical patterns are shown in the inset plots. We measure 
the probability of propagation by determining whether infected 
individuals remain after five hundred steps. By averaging over one 
hundred runs, we obtain a phase transition diagram displaying a 
sharp change. For our example, where the probability of recovery 
is fixed to 1/2, the simulation reveals that epidemic spreading 
does not occur if the rate of infection is smaller than 1/2 whereas 
a stable infected population is always observed otherwise. This 
implies the presence of a second-order (critical) phase transition 
and the existence of a threshold of infection rate that needs to 
be overcome in order to generate a stable infected population.” 
A percolation process is at work. 


' Alternatively, the process can also model the occupation of a given space by 
a population under a stochostic process involving. 

? This is actually an example of the absorbing phase transition mentioned at 
the end of chapter 1. The absorbing phase corresponds to the zero infection 
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Figure 9.1. Epidemic spreading in a model based on a contact 
process. Here a one-dimensional lattice is considered, with each site 
occupied by an infective (gray ball) or a susceptible (empty ball). In 
(a) we show the four possible transitions that can take place: either 
an individual becomes susceptible (independent of the state of its 
neighbors) or a susceptible becomes infected. Three possible situa- 
tions are allowed. If we measure the fraction of infectives after a long 
transient for different infection levels (b), the resulting curve displays 
a transition from no infectives to a sustained epidemic state (see text). 


A theoretical treatment of the process including the spatial 
correlations is difficult (Hinrichsen 2000; Marro and Dickman 
1997). If p(x, £) indicates the probability of infection at a given 
point x € £, that is, 

p(x, t) = P[S,(¢) = 1] (9.1) 
then p(x, ¢) will evolve (according to the previous rules) follow- 
ing: 

q 
a = —pp(x, t) + ei py P [S(t) = 0, Sa(t) = 1] (9.2) 


<u> 


phase. Once reached, there is no escape toward the second phase. Of course this 
is true under the hypothesis that the system is isolated and thus no infections can 
come from the external environment (which is clearly not the case in general). 
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where <u> indicates the sum over the set of g nearest neighbors. 
For one dimension, this set is simply {2 — 1,2 + 1}. The 
previous equation is exact, but its computation would require 
knowledge of the probabilities associated with the interactions 
between nearest sites. As will be shown below, a qualitative 
understanding of the problem is easily achieved from a simple 
mean field argument. 


9.2 SIS Model 


Here we define a minimal model of epidemic spreading, taking 
place between two populations, as previously defined for the 
contact process (see Murray 1990 and references therein). For 
simplicity, we will also indicate their population sizes as S and J. 

The infection dynamics will include two basic processes: (a) 
new infections occur provided that infected individuals contact 
healthy ones. The rate at which infection takes place will be 
indicated as u. On the other hand, (b) infected individuals can 
recover at a rate œ, becoming healthy again. The two components 
of the process can be summarized in terms of two chemical-like 


reactions: 
Pep 3 (9.3) 
for the infection event and 
(25% (9.4) 


for the recovery process. Additionally, let us assume that the 
total population size N remains constant, that is, S + J = N. 
As a result, we can write down the equation for the number of 
infectives as 


dI 
gy TYS = N-I) -oal (9.5) 


As we can see, we have two conflicting components defining the 
infection dynamics. One is a density-dependent term where the 
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Figure 9.2. Epidemic spreading in the real world. In (a) we show 
the number of months with measles as a function of population 
size. Note the presence of a sharp change at intermediate population 
values. A simple model of epidemic spreading (b) considers two sub- 
populations composed of infected and healthy individuals. Although 
infected ones recover with some rate œ, newly infected ones require 
interaction between both groups. 


two groups need to interact, whereas the second is the linear 
decay of infected individuals, as they recover from infection. This 
equation has two equilibrium points, namely /* = 0, describing 
the all-healthy population state and a second one 


on al (9.6) 
H 
describing a population with a finite number of infected individ- 
uals. It is easy to see that the stability of these states is determined 
by the presence of a critical condition defined by the so-called 
basic reproductive rate of the pathogen: 


Ppa ad (9.7) 
Q 


Epidemic spreading will occur if Rọ > 1, and decay will take 
place instead for Ro < 1. It is interesting to note that Rọ involves 
several components, including the infectivity of the pathogen ju 
but also the population size NV. This implies that as population 
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Figure 9.3. Time series of epidemic spreading. Here we use a = 1.0 
and (a) u=0.8 (subcritical) and (b) u = 1.2 (supercritical). The 
inset in figure (a) is a linear-log display of the trajectories at the 
subcritical phase, showing a linear decay that corresponds to an 
exponential behavior. 


size increases so do the chances for a pathogen to spread. This is 
actually a well-known observation, and an example is displayed 
in figure 9.2a for measles. Here we plot a surrogate of the 
previous order parameter using the (average) lifespan of measles 
outbreaks in different populations. Below some threshold, no 
effective spreading occurs, whereas after M ~ 30000, spreading 
seems guaranteed. 


In order further analyze our model, we will use a normalization of 
population sizes, that is, x = Z/N and y = S/N and thus using 
y = 1l — x, the previous equation reduces to: 

dx 


a px(1 — x) —ax (9.8) 


and the equilibrium points now read xj = 0 and 


geis 5 (9.9) 
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Figure 9.4. Phase transitions in the epidemic spreading model. In 
(a) the bifurcation diagram for the mean field model is shown, using 
a = 1.0, giving u, = 1.0. In (b) the two phases in parameter space 
(a, u) are shown. 


respectively. The epidemic will die out if 
Ay(0) =w-—a <0 (9.10) 


that is, if > u. The opposite case, defined for œ < u corresponds 
to a stable epidemic spreading. The two points exchange their 
stability. In summary, as shown by the bifurcation diagram in 
figure 9.4a, the critical point ue = @ separates the two phases of 
this model: (1) the subcritical phase, where the epidemics dies 
out and (2) the supercritical phase, where the epidemics self- 
maintains. The corresponding parameter space is also displayed 
in figure 9.4b. The results essentially tell us that the balance 
between infection and recovery define the limits separating the 
two possible qualitative outcomes of the model. 


9.3 Vaccination Thresholds 


The presence of a critical threshold in epidemic dynamics implies 
that tipping points exist. Once a pathogen is able to cross the 
boundary, epidemic spreading will occur. But contention based 
on vaccination strategies can also exploit this threshold behavior 
in efficient ways. Consider our population, but now imagine that 
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Figure 9.5. Effects of vaccination on epidemic spreading. Here 
different levels of vaccinated individuals R are introduced and 
compared with the predicted curve for nonvaccinated population 
(continuous line). 


a fraction R of its individuals has been vaccinated and thus cannot 
be infected. Assuming normalized population size (NV = 1), the 
new equation will read: 

dx 


— = ux(l — R — x)— &æx (9.11) 
dt 


with 0 < R < 1. The reader can easily check that the new steady 
states for the infection phase are now 
a 
alka (9.12) 
u 


which are stable if 


u> WR) = — (9.13) 


The resulting bifurcation curves are displayed in figure 9.5, 
where we can see that the eradication threshold shifts to higher 
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values as the R parameter is increased. An alternate form for this 
threshold (as the reader can check) can be written as: 


1 

R>1- R (9.14) 
(see Anderson and May 1991). The target of any eradication 
policy requires reducing Ro to below one. Available data allows es- 
timating this parameter. For smallpox, for example, we have Ro © 
3 — 5, and 60-70 percent of the population must be vaccinated. 
However, measles has a much higher basic reproductive rate of 
Ro œ% 13, which requires a much larger vaccinated population 
(around 20 percent or R ~ 0.92). It is worth noting that smallpox 
has been eradicated from the world, and that measles would be 
much harder to eradicate. 


9.4 Pathogens and Networks 


The SIS model approach is to course a first approximation of 
epidemics, although it allows us to understand many interesing 
cases. But the population biology of infectious diseases involves 
further complexities (Anderson and May 1991). These include 
spatial dynamics (Grenfell et al. 2001), evolutionary change 
(Earn et al. 2002), and the structure of human-based networks 
(Pastor-Satorras and Vespignani 2000; May and Lloyd 2001; 
Colizza et al. 2007). As mentioned at the end of chapter 5, 
the organization of real networks is usually heterogeneous. This 
includes several related webs, particularly two: transportation 
networks and social interaction networks. The first group is 
dominated by airports: a given individual can travel from one 
part of the globe to another very easily. Most airports have one 
or two links, and a few dominate the scene. These hubs connect 
to many parts of the world and they largely influence the like- 
lihood of a given infectious disease to propagate worldwide as a 
pandemic. Moreover, on the social level, sexual contact networks 
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have been shown to display scale-free organization too (Liljeros 
et al. 2001). 

The implications of this pattern have been widely discussed 
because of their potential use in controling the spread of both 
biological and computer viruses (Lloyd and May 2001). In a 
nutshell, both large airports and promiscuous sexual partners need 
to be taken into account. The first is crucial because it allows 
pathogens to jump to multiple, geographically distant targets. The 
second is significant because it represents the highest risk group, 
given its leading role in allowing pathogens to be transmitted. The 
presence of this risk group modifies one of our main findings: it 
reduces or even removes the presence of an eradication threshold 
(Pastor-Satorras and Vespignani 2001). This lack of a phase 
transition implies that epidemics will always propagate due to the 
many opportunities provided by highly connected nodes. We can 
return again to our original scenario with a threshold by properly 
influencing the hubs (Dezso and Barabdsi 2002; Pastor-Satorras 
and Vespignani 2002; Zanette and Kuperman 2002; Latora et al. 
2006; Colizza et al. 2007). 
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GENE NETWORKS 


10.1 Genes and Cell Types 


Any multicellular organism, from sponges to humans, is formed 
by a vast amount of cells defining an integrated being (Wolpert 
2004). Cells are the fundamental units of life, and within a com- 
plex organism, they cooperate in different ways, forming tissues 
and organs and sustaining the life of the individual in a reliable 
and homeostatic manner. From a single initial cell, a multicellular 
organism develops into a complex life form with different cell 
types (figure 10.1a). Cell types emerge through the process of 
differentiation (Alberts et al. 2004). Such process involves the 
interactions between a cell’s gene regulatory networks and the 
information sent by neighboring cells and the environment. 
How are different cell types generated? Mounting evidence 
points toward genetic switches of different types as the underlying 
mechanism of cell differentiation (Huang et al. 2005). Switching 
mechanisms would be responsible for creating distinct patterns 
of gene activation and repression. If cells are in physical contact 
or can communicate through diffusion, nonlinear mechanisms 
of cell-cell information exchange can propagate signals through 
space and lead to a stable spatial pattern. The regulatory pattern 
of interactions is expected to be complex and involve both positive 


Figure 10.1. Tissues and organs are formed by combinations 
of different cell types (a) displaying very diverse morphologies. 
Each type of cell displays a number of specific functionalities 

(continued) 
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and negative links. An example of a given subnetwork is shown in 
figure 10.1b (Espinosa-Soto et al. 2004). This network involves 
regulatory interactions associated with the determination of cell 
fate in the flowering of a species of plant, Arabidopsis thaliana. 
Untangling the wires of such complex networks is one of the main 
challenges in systems biology (Bornholdt 2002, 2005; Alvarez- 
Buylla et al. 2007). 

This problem involves many degrees of freedom. Moreover, 
it is remarkable that the possible and the actual in the universe 
of potential cell types are very different. Each cell type can be 
understood (on a first approximation) as a set of ON-OFF 
genes: some genes will be activated, generating special proteins 
characteristic of the given cell type, whereas others will be 
inactivated (figure 10.1b). Each cell shares the same set of NV 
genes (from thousands to tens of thousands) and thus the potential 
number of cell types N; under this ON-OFF picture will be 


astronomic, namely: 
N, ~ 2% (10.1) 


whereas in reality it remains confined to a much smaller number, 
from tens to hundreds (Kauffman 1993; Carroll 2001). Why 
is this number so small? More important, cell types are stable 
states, not easily able to switch from one to another without 


Figure 10.1 (continued). and cooperation among cells allow tissues 
and organs to perform optimaly. Each one seen, as a first approxima- 
tion, as the result of the activation of a given set of genes (drawing 
by R. Solé). In (b) we show an example of a gene network involved 
in the determination of cell states of the floral organ of the plant 
Arabidopsis thaliana. Here network nodes represent active proteins 
of corresponding genes, and the edges represent the regulatory 
interactions between node pairs (arrows are positive, and blunt- 
end lines are negative). Dashed lines are hypothetical interactions. 
(Redrawn from Espinosa-Soto et al. 2004) 
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some perturbation. On the contrary, cellular states tend to be 
stable over time, consistent with the desired stability for a complex 
organism to survive. 


10.2 Boolean Dynamics 


How to deal with the dynamics of cells at the genetic level, if the 
number of genes and attractors is so high? One possible shortcut 
is given by Random Boolean Networks (RBN), a class-of-discrete 
dynamical system where each element is synchronously updated 
through a specified function: 


Si + 1) = Pls), Sa), ox dye) (10.2) 


where ¿ = 1,..., Nand S;(¢) € & = {0, 1} are the only two 
possible states (ON-OFF). Let us note that each element has a 
different associated function, since different sets of genes send 
input to it (that’s why we use an index for each ® function). 
The previous definition assumes that time is discrete (states are 
updated using some clock having a characteristic time scale) and 
that each gene receives input from (is regulated by) exactly K 
other genes (where K = 0,1,...). Many potential choices 
can be made in defining the exact Boolean functions ®;. The 
easiest approach is to choose each randomly from the set of all 
K-inputs Boolean functions. In other words, for each K-string 
of possible inputs, we randomly choose the output (0 or 1) with 
some probability p. The links between elements are also chosen 
at random (as in a random graph, chapter 5). 

Gene networks able to generate multiple attractors have been 
identified in real systems (Buylla et al. 2007). An example is given 
in figure 10.1b, where we display the gene-gene interaction map 
associated to early flower development (Chaos et al. 2006). In this 
case, six cell types are matched by six attractors obtained from 
a Boolean representation of gene states. These models exhibit a 
range of dynamical patterns, from fixed points to oscillations and 
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chaos (Aldana et al. 2003; Drossel 2007). Fixed points are sets 
S* = (ST, ..., Sy) such that 


SEE] (10.3) 


for all ¿ = 1,..., N. A simple example of these networks is 
shown in figure 10.2a-b, where a small N = 3, K = 2 case 
is being considered. Here the functions used are AND and OR. 
These Boolean functions are simply defined by means of Boolean 
tables, namely: 


So(t) | S) | SE 1) S) | S) | So(¢+ 1) 
0 0 0 0 0 0 
0 1 0 0 1 1 
1 0 0 1 0 1 
1 1 1 1 1 1 


for the AND and the OR gates associated to Sı and S2, respec- 
tively. 

Using these tables, it is easy to compute the next state exhibited 
by the whole network (Kauffman 1993). Even in this simple 
example, we can see that the possible dynamical patterns are 
nontrivial. In figure 10.2c we show the possible transitions 
between different states. If the system starts at (0, 0, 0), it remains 
there (point attractor) whereas if it starts from (0,0, 1) or 
(0, 1,0), these two states flip into each other defining a cycle. 
The rest of the states flow toward the second fixed point (1, 1, 1) 
thus defining a basin of attraction. 


10.3 Percolation Analysis 


The previous example is just a picture illustrating how the 
dynamics is defined. When N becomes large, the combinatorial 
possibilities rapidly explode and it does not make much sense 
to use tables. However, from a statistical point of view, two 


Figure 10.2. Boolean networks and attractors. In (a) a hypothetical 
gene network is shown, where each gene (long rectangles) produces 
a protein (gray balls) that regulates the other two genes. An abstract 
formulation of this small genetic circuit is shown in (b) where the 
explicit mechanism of gene-gene interaction is replaced by arrows 
indicating the existence of an interaction (here N = 3 and K = 2 
is used). For a given choice of the Boolean functions (an OR gate 
for genes 2 and 3 and an AND gate for 1, see text), the possible 
transitions between states for the state space {5S}, S2, S3} are shown 
in (c). Three attractors are obtained. 


10 


order 


Figure 10.3. Phase transition in random Boolean networks. The top 
left plot displays the two phases in the connectivity-bias parameter 
space. Three examples of the spatiotemporal behavior are indicated 
in (a—c) for three different values of K at p = 0.5. Here time goes 
from left to right, and black and white squares indicate 1 and 0 states, 
respectively. 
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well-defined phases have been identified (Derrida and Stauffer 
1986; Derrida and Pomeau 1986; Kurten 1988; Kauffman 1993; 
Glass and Hill 1998; Aldana et al. 2003; Drossel 2005). Assuming 
that the probability of having a zero or one in the Boolean table 
is the same (i.e., we generate the two states using a fair coin), 
we can tune the number of connections K. For K =1, the 
system is frozen: every element rapidly achieves a given steady 
state and remains there. For K > 3 the system is chaotic, with 
elements changing in a seemingly random fashion. This can be 
further generalized by introducing an additional parameter p that 
measures the probability that the output of equation (2) is zero. 

These two phases describe two different scenarios involving 
order and disorder. Moreover, the attractors in each phase are 
highly sensitive to perturbations: in both cases, a change in the 
state of a given element generates changes in a system’s state. 
At the ordered domain, an element that is randomly chosen and 
changed remains in the same state (the input from others is too 
weak and unable to reverse the flip). Instead, in the disordered 
phase, a flip triggers a chain reaction that propagates through the 
entire system (too many connections). In both cases, a very large 
number of attractors is present. However, for K = 2 the system 
displays a reduced number of attractors that is highly stable: a 
given perturbation is typically reversed and thus homeostasis is at 
work. Such properties are compatible with what we would expect 
for a cell type, as discussed earlier. 

The dynamical equations describing the RBN model (2) can 
be seen as a set of tables. For each input, we have a given output: 


bn (2), Sel Dy aces hep) — > Se +e 1) (10.4) 


and each value S;(¢ + 1) (i.e., the right column in the table) is 
either zero or one with some probability p. In other words, we 
have, for any S; in the system, 


P[®;(S;, (t), Sa (2), EERS S;,(t)) = 1] = p (10.5) 
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and 
PEP; (SiC), Salt), n Si) = 0) =1= p (10.6) 


This parameter is the so-called bias and it seems clear that it 
will also influence the location of the phase transition.! In our 
previous example, we used p = 1/2, but a very small (or large) p 
will clearly make the system much more homogeneous and thus 
more likely to be ordered, no matter the value of K. 

In order to derive the mean field phase transition line, we 
need only consider a simple percolation argument (Luque and 
Solé 1997). At the subcritical phase, a change in a given element 
is unable to propagate through the system, whereas within the 
chaotic phase any change has an impact. This reminds us of the 
percolation phenomenon discussed in chapter 4: provided that 
a critical threshold is reached, a single flip in the network can 
generate a cascade of changes. In our context, a flip means that 
we change one input element, namely: 


(S, wz AE) vt eg Si) e Sa Oo vee LG). 2x5 Si 
(10.7) 


What is the probability that such change generates further 
changes in the other elements that receive input from S;? In figure 
10.4a the basic situation is shown using a simple tree structure. 
Here at the bottom we have the element that is in a given state 
(1 in our example). This element send inputs to others (an average 
of (K)) which are also in a given state. These send further inputs 
to others. If we flip this element (i.e., 1 — 0, figure 10.4b), what 
is the likelihood that this will trigger an avalanche? 

If we choose an arbitrary element, the question is what is the 
probability that the output in the table defined by ®; will change 


! For the small examples defined above and used to generate the attractors of 
figure 10.2, we have p(AND) = 1/4 and p(OR) = 3/4. But no extrapolations 
from these small-sized examples can be made in terms of the global behavior of 
large systems. 
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Figure 10.4. Propagation of perturbations in random Boolean net- 
works. In (a) we show part of a given network using a treelike 
representation, downstream from a given element. This element 
(bottom) send inputs to K other elements which send inputs to 
others. If we flip the state of the bottom element (b) further 
changes can be induced, provided that the Boolean tables allow 
this to happen. In (c) we plot in black the propagating changes, 
which can percolate if the appropriate critical condition is reached 
(see text). 


to a different value. This is easily found: if the initial output is 1, 
the probability that it changes to zero is simply p(1 — p). Since 
we need to consider the symmetric situation (0 — 1), we actually 
have a total probability of change 2 p(1— p). Every element affects 
(K) others, and thus the total number of changes reads 


N = 2p(1— p)(K) (10.8) 


Now we are ready to find the critical line: the boundary separating 
propagation from a frozen state is simply N = 1, which gives: 


1 
K: =— 10.9 
= yaa (10.9 
and for the special case p = 1/2, we have K,(1/2) = 2. 
The previous formalism can be easily applied to more complex 
situations, including random networks with multiple states and 
neural networks (Luque and Solé 1997). 


118 Gene Networks 
10.4 Is Cell Dynamics Boolean? 


The assumption that gene network dynamics can be represented 
in terms of Boolean, discrete sets of states, is of course a strong 
one. Such an approach is necessarily a crude approximation of 
real genes. But we also need to keep in mind that the level of 
approximation used is directly related to the questions we are 
trying to answer. 

The correctness of the ON-OFF picture has received con- 
siderable attention over the years (see Bornholdt 2008 and 
references cited; Macia et al. 2009). The truth is that the logic 
of gene regulation has much in common with a computational 
device (Bray 1995; Hogeweg 2002; Macia and Solé 2008), and 
switchlike genetic behavior is known to be more the rule than 
the exception. In order to illustrate this point, let us consider the 
simplest example of a genetic switch, using a single gene coding 
for a protein P that activates the gene through a positive feedback 
loop (figure 10.5a). Such a scenario has been shown to occur 
experimentally by tuning the strength of transcription of a given 
gene (Isaacs et al. 2003). The protein is assumed to form a dimer 
in order to perform its regulatory function. If p indicates the 
concentration of P, it can be shown that the kinetic description 
of this system is given by the nonlinear equation: 


dp _ ap? 


dt 1+2? 


— 3p (10.10) 


where ô is the degradation rate of monomers. 
Together with a trivial fixed point p* = 0, which is always 
stable, two additional points are obtained, namely 


aaa G -4 (10.11) 


The higher is stable whereas the lowest is unstable. If we use 
&œ/ô as a bifurcation parameter, we can see that these two fixed 
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Figure 10.5. In (a) a simple example of a genetic circuit displaying 
two possible alternative states is shown. Here a gene is expressed and 
translated into a protein P (gray balls). This protein degrades and is 
removed (empty ball) or forms dimers which can bind to the DNA 
and trigger further synthesis. Experimental results using engineered 
systems (b) confirm that this system is able to behave as a switch 
(adapted from Isaacs et al. 2003). 


points will exist provided that œ > 26. In other words, two basic 
alternative states exist if the previous inequality holds. The one 
chosen might depend on how other cellular signals tune @ or ô. 

The previous example illustrates our point: although strictly 
speaking the underlying system is continuous, the nature of the 
final decisions can be understood in terms of switching logic. 
Moreover, we might ask why bistable, Boolean dynamics is so 
widespread: why not tri-, four-, or n-stable? Although tristable 
dynamics might be involved in some special circuits, such as 
stem-cell switches (Huang et al. 2007), it seems that, beyond two 
alternative states, switches become less likely to properly work in 
a reliable way (Macia et al. 2009). We could say that nature has 
no choice but being (nearly) Boolean. 


11 


CANCER DYNAMICS 


11.1 Abnormal Growth 


Cancer is the result of a system breakdown that arises in a cell 
society when a single cell (due to a mutation or set of mutations) 
starts to display uncontrolled growth (Hanahan and Weinberg 
2000; Weinberg 2007). The cooperation that maintains the 
integrity of a multicellular organism is thus disrupted. Further 
changes in the population generated by such abnormal cell growth 
can lead to malignant tumor growth, eventually killing the host. 
From an evolutionary point of view, tumor progression is a 
microevolution process in which tumors must overcome selection 
barriers imposed by the organism (Merlo et al. 2006). The 
emergence and evolution of tumors involve a number of phe- 
nomena that are well known in complex systems (Kitano 2004; 
Anderson and Quaranta 2008). These include, for example, 
growth, competition (Michelson 1987; Gatenby 1996), spatial 
dynamics (González-Garcia et al. 2002; Cristini et al. 2005), and 
phase transitions (Garay and Lefever 1977; Solé 2003; Solé and 
Deisboeck 2004). 

As discussed in (Alberts et al. 2007), a multicellular system is 
not far from an ecosystem whose individual members are cells, 
reproduced in a collaborative way and organized into tissues. In 
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this sense, understanding multicellularity requires concepts that 
are well known in population dynamics such as birth, death, 
habitat, and the maintenance of population size (Merlo et al. 
2006). Under normal conditions, there is no need to worry 
about selection and mutation: as opposed to the survival of the 
fittest, the cell society involves cooperation and, when needed, 
the death of its faulty individual units. Mutations occur all the 
time but sophisticated mechanisms are employed in detecting 
them and either repairing the damage or triggering the death of 
the cell displaying mutations (Weinberg 2007). Abnormal cells 
can be indentified from within (i.e., through molecular signaling 
mechanisms operating inside the damaged cell) or by means 
of interactions with other cells. The latter mechanism involves 
immune responses. 

Selection barriers (such as an attack from the immune system 
or physical barriers of different types) can be overcome by a tumor 
provided that the diversity of mutant cells is high enough to 
generate a successful strain. High mutation rates are thus a way to 
escape from the host responses and it is actually known that most 
human cancers are genetically unstable (Weinberg 2007). Genetic 
instability results from mutations in genes that are implicated in 
DNA repair or in maintaining the integrity of chromosomes. As 
a result, mutations accumulate at very high rates. 

Modeling cancer is far from trivial (Wodarz and Komarova 
2005; Spencer et al. 2006; Byrne et al. 2006; Anderson and 
Quaranta 2008). Although the first impression is that tumors are 
simply growing, disordered pieces of tissue, the reality is that they 
involve several levels of complexity. Cancer is a complex adaptive 
system and as such exhibits robustness (Kitano 2004; Deisboeck 
and Couzin 2009) and thus adaptability to changing conditions. 
However, as any other complex system, cancer can also display 
Achilles’ heels. In this context, models offer the possibility of 
exploring potential sources of fragility such as catastrophes and 
breakpoints. Two examples are presented below. 
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11.2 Tumor Decay under Immune Attack 


As mentioned above, cancer cells must overcome different types 
of selection barriers in order for tumors to expand. Among these 
barriers, the immune system (IS) appears as a specially important 
one. The immune system is itself a complex system, able to 
respond to the invasion of foreign pathogens with high efficiency 
(Janeway et al. 2001). There is no space to provide even a basic 
overview of how the IS works. We will confine ourselves to a 
rather important component and how it interacts with cancer. 
In a nutshell, the IS response is based on the recognition of a 
very large class of molecules, also known as antigens. These are 
mostly carried by infectious agents, such as viruses and bacteria.! 
Immune responses are directed toward detecting and eliminating 
these molecules and their carriers. However, if the carrier is able 
to change itself (such as RNA viruses, see chapter 6), it can escape 
from the immune attack. Different escape strategies are known to 
occur, all requiring some degree of internal plasticity. 

Among the key mechanisms at play, an important part of 
the IS response is based on special types of cells able to detect 
anomalous molecular markers at the surface of infected cells. 
These include the cytotoxic T cells (CTLs). Each class of CTL 
is able to recognize one specific antigen among all the huge 
set of possible molecules. Once recognition occurs, an attack is 
launched and the CTL injects toxic granules into the target cell 
(figure 11.1) eventually killing it. 

In principle, these mechanisms could be a source of cancer- 
cell removal. However, tumor cells can evade IS recognition in 
different ways (Weinberg 2007), and sometimes the IS becomes 


! However, the IS needs also to distinguish between nonself and self: it must 
avoid attacking normal cells. The IS actually learns to distinguish what is 
part of the body, but such recognition is not perfect and mistakes can ocur. 
Such mistakes happen from time to time causing different types of so-called 
autoimmune diseases. 
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Figure 11.1. Attack on a cancer cell unleashed by a CTL lympho- 
cyte, as shown in (a) where the larger cell belongs to a neoplasic tissue 
and the small spheres are cytotoxic cells. In (b) the sequence shows 
the whole process, including the formation of a cell-cell complex and 
the final killing of the tumor cell (adapted from Weinberg 2007). 


ineffective against cancer growth. The interaction between cancer 
and the immune system actually involves multiple layers of 
complexity (Adam and Bellomo 1997). An especially important 
problem is determining the conditions under which the IS 
will be able to remove cancer, and some models reveal the 
presence of transition phenomena. Below we present the Garay- 
Lefever model of cancer—immune system interactions (1977), 
which is shown to display both first- and second-order transition 
phenomena (Garay and Lefever 1977). 

Consider a population of tumor cells, to be indicated as X, that 
can be attacked by some class of cytotoxic killer cells (T-cells). 
Here cancer cells replicate at a given rate 7. In this model, T cells 
can recognize and kill cancer cells by detecting special markers on 
the tumor cell surface. The identification occurs at a rate $; and 
triggers the formation of a cell-cell complex C (figure 11.1). Once 
the cancer cell has been killed, the complex dissociates and the T 


cells are free again. 
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The previous components define a set of reactions, namely: 
X> 2X (11.1) 
for cancer growth kinetics and 
Tex“ C4 yao (11.2) 
for cancer-CTL interaction. Here D indicates a dead (cancer) 
cell. The rates k; and k measure the speed at which the 


complex forms and disintegrates, respectively. The corresponding 
equations describing this dynamics are: 


ge Xil a ky TX (11.3) 
a K : 
T 
dt 


where a maximum number of cancer cells K is assumed to 
exist. The first term on the right-hand side of the cancer- 
growth equation is thus consistent with a logistic growth model 
(chapter 2)? 

As defined here, cytotoxic cells act as predators, reducing the 
number of cancer cells at a rate kı TX, that is, depending on the 
encounter rates, efficiency of recognition, and other phenomena 
captured by kı. Moreover, killer cells forming complexes get back 
into the free T cell pool as soon as the cancer cell is killed at a 
rate k>. 

This two-dimensional system can be reduced by further assum- 
ing (Garay and Lefever 1977) that the lifetime of the C complex 
is short and thus a steady state assumption: 


— 0 (11.5) 
can be made. This equilibrium condition just tells us that 


cytotoxic cells form complexes that rapidly split (compared to 


? Actually, a more realistic representation of tumor growth dynamics is 
provided by the so-called Gompertz model, where the kinetic equation is given 


by dx/dt = nX log(K/X). 
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other phenomena involved in the model). The total number of 
cytotoxic cells, £, can be computed as follows 


E=T+C (11.6) 


since E includes both complex-associated and free killer cells. In 

the following, we assume that Æ is a conserved (constant) value. 
Using the previous stationary approach for T cell dynamics, 

we have 7(X) = k2C/kıX and using C = E — T, we obtain: 


T(X) = (11.7) 


1+4xX 


which allows us to write down the whole mean field equation for 


tumor growth dynamics in terms of a one-dimensional system,’ 
namely: 
dX X kıX E 
uta d G Ea (11.8) 
dt K 1+ Ex 


Following Garay and Lefever, we will reduce the number of pa- 
rameters in this model by using the following rescaled quantities: 
u = ki E/n, 0 = ko/ki Nand x = k,X/k (time t would be 
rescaled accordingly as t + nt, but we keep the same notation). 
The model now reads: 

dx [Lx 


ee ee (11.9) 


Using this reduced one-dimensional model, three possible fixed 
points are obtained (figure 11.2), namely x* =0 (cancer extinc- 
tion) and the two nontrivial solutions: 


1 
= = = . 
a= E 0/0 +0} 40u] (11.10) 


The tumor-free fixed point will be stable provided that 


Ox 


àn (0) = (=) <0 (11.11) 
x=0 


3 The resulting model is very similar to the one obtained by considering the 
interaction between cancer cells and a given anticancer drug. 
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Figure 11.2. Bifurcations in the Garay-Lefever model. Here we plot 
the stable solutions (thick, continuous lines) for (a) 0 = 1.5 and 
(b) 0 = 0.5 corresponding to the monostable and bistable domains 
(see text). 


which in our case is simply u > ue = 1. Since the parameter u 
captures the conflict between the rates n (tumor cell formation) 
and kı (tumor cell removal) the critical condition w,=1 natu- 
rally separates cancer survival from extinction. 

Looking at the nontrivial solutions, we can see that two 

possible scenarios can be defined: 

1. For 0 > 1, the fixed point x_ will be negative (and thus 
biologically irrelevant) whereas x4 remains positive 
provided that u < 1. In this case, a cancer cell 
population will exist below the critical point pe, 
exchanging stability with the trivial solution. 

2. For @ < 1, when the inequality (1 + 6)? > 40(u — 1) 
holds, the two nontrivial solutions exist. It is 
not difficult to see that, given the previous results, the 
interval of existence of these two simultaneous 


solutions is: 
i (+8) 
zye ——— 
i 40 


defining a domain of bistable behavior. For u < 1 only 


=i (11.12) 


the first solution x+ is present. For the bistable domain, 
we can describe the two solutions as follows: 
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Figure 11.3. The three phases exhibited by the Garay-Lefever model 
(1977). Here the upper white domain indicates that no cancer cells 
are stable (tumor extinction) whereas two other areas correspond 
to the presence of monostable cancer (white low domain) and 
two alternative states (gray area). Four examples of the associated 
potentials are also shown, with filled circles indicating the presence 
of stable states. 


1 
mas [1 ~6+ fon »)| (11.13) 


In figure 11.2 we show two examples of the bifurcation diagrams 
associated with the model. In (a) an example of the continuous 
transition scenario is obtained for 89 > 1 whereas the bistable 
case, for 0 < 1, is shown in (b). In figure 11.3 the basic results 
are summarized by representing the three phases on the (0, p) 
parameter space. We also plot four examples of the associated 
potential function, here given by 
2 3 


x Ox 
®,,(x) = -F + Ey + p(x — ln(1 + x)) (11.14) 


which also illustrate the exchange of stability between different 
attractors. 

One especially interesting result of this model is the existence 
of a catastrophic shift from cancer to cancer-free tissue. Looking 
at the transition shown in figure 11.2b, we conclude (under our 
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special assumptions) that a slight increase in jz could trigger 
the remission of cancer and eventually its disappearance. Such 
phenomena have actually been reported in some rare cases where 
advanced (usually lethal) tumors were shown to experience spon- 
taneous reversal (Cole 1981). These spontaneous remissions have 
also been observed in some animal models (Cui et al. 2003). 
The mechanisms pervading the dynamics of these remarkable 
events are largely unknown, but a late, powerful immune response 
is likely to be involved (see Cui et al. 2003; Weinberg 2007). 
The almost sudden regression experienced by these subjects 
seems consistent with the first-order shifts displayed by our toy 
model. 


11.3 Thresholds in Cancer Instability 


In this section we consider a different scenario in which thresh- 
olds to cancer viability might exist. The key property now is 
genetic instability. More precisely, we consider the growth of 
tumor populations in those types of cancer where mechanisms 
involved in the maintenance of genomic integrity have failed 
(Lengauer et al. 1998; Hoeijmaker 2001; Weinberg 2007). Such 
failure is linked to mutations in a class of genes responsible for 
detecting and controlling potential mistakes in the process of 
DNA replication. This leads to the so-called mutator phenotype 
(Loeb et al. 2003). If DNA stability is jeopardized, increasing 
levels of instability foster tumor progression and further instability 
(see Cahill et al. 1999; Loeb et al. 2003). This leads to the 
question of what levels of instability can be achieved before 
harmful mutations stop further progression. Such instabilities are 
observable at different levels, from chromosomal aberrations to 
pathological tissue organization. The latter is well illustrated in 
figure 11.4, where different stages of colon cancer are shown 
at different times of tumor development (Gatenby and Frieden 
2002). A clear loss of order can be appreciated, although some 
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Figure 11.4. Degradation of tissue organization at different stages of 
the evolution of a colon cancer (adapted from Gatenby and Frieden 
2002). Here we have: (a) normal tissue (b) colon adenomatous 
polyps with low grade, (c) dysplasia, and (d) undifferentiated, 
invasive colon cancer with highly disordered morphology and loss 
of normal colon tissue organization. 


degree of coherence is observable as cells maintain some capacity 
to differentiate. Eventually (as illustrated by the last picture in 
the sequence) a complete loss of structure occurs. It is known 
that many of these cancer cells are actually not able to repli- 
cate, thus suggesting that (as it occurs with RNA viruses, see 
chapter 7) some thresholds of viability might be present (Cahill 
et al. 1999; Solé 2003; Solé and Deisboeck 2004). If true, then 
such thresholds could be exploited for cancer treatment. Below 
we present a minimal model of cancer growth illustrating these 
ideas. 

Assuming that two populations of cells are at play, namely 
normal (7) and cancerous (x), we will introduce a simple, two- 
dimensional model of cancer instability. The model assumes that 
both populations are homogeneous, thus all cells in each subset 
are identical in their replication capacities. This is of course 
an oversimplification, particularly when dealing with unstable 
cancer, where by definition high diversity is expected. Assuming 
that these populations compete for available space and resources, 
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the simplest model reads: 


ee (11.15) 
dt 
a eee oe (11.16) 
dt 


where 7 is the replication rate of normal cells in the original 
tissue and ju indicates the rate of genetic instability.4 The function 
T (u) includes the (nonlinear) effect of instability on the growth 
of unstable cancer cells. Although u can be understood as a 
mutation rate, other changes beyond point mutations are being 
taken into account (such as chromosomal rearrangements). The 
last terms on the right-hand side of both equations introduce an 
outflow that is required to guarantee a finite population size. 

For simplicity, we assume that there is a maximum population 
size and that the total population is constant, that is, n + x = C, 
which for simplicity we fix to C=1. This is the constant 
population constraint used in chapter 6, which implies that 
dx/dt + dn/dt = 0. Using this constraint, we obtain: 


F(n,x) =rn+rT(u)x (11.17) 
And since n = 1 — x, we can reduce our previous model to the 
following one-dimensional system: 


n rT (u) — 1)x(1 — x) (11.18) 
dt 


which is nothing but a logistic equation (chapter 2). We know 
from our previous analysis that two fixed points are present: the 
zero-population one x* = 0 and the maximum population state, 
here x* = 1. It is easy to see that the first is stable if r (u) < 1 
and unstable otherwise. By properly defining the function T (u) 
we might be able to define the conditions under which genetic 


í Previous models used a formal approach based on the quasispecies model 
developed in chapter for RNA viruses. 
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instability allows cancer growth to occur and overcompete the 
host tissue. 

In order to propose a reasonable form for T (u), two features 
of instability must be considered. The first is that, at low levels, 
instability allows the emergence of clones able to grow faster and 
overcome selection barriers. In other words, some amount of 
instability is good for the cancer population: it provides a source 
for adaptability and selective advantages. The second deals with 
the negative part of the story: at high levels, mutations tend to 
be harmful and cells nonviable. The function should be such that 
r(0) = 1 (i.e., no effect is at work) whereas we should expect 
T(t) —> 0 as instability grows (i.e., u > 00). One possible choice 
(Solé et al. 2008) is using: 


Tu) = (1 + aue À (11.19) 


which is consistent with all the previous requirements. Here œ > 0 
is a parameter that weights the selective advantage of genetic 
instability: the higher its value, the faster the increase in growth 
rate. The parameter u* is simply some characteristic rate that 
we fix here to u* =0.5. By using this functional form for T (u) 
we find two phases, as indicated in figure 11.5. Since the two 
alternative possibilities (x* =0 and x* = 1) are separated by a 
critical line, a sharp change is expected to take place near this 
boundary. The corresponding potential is now: 


2 3 
(x) = -rT (u) — 1) (= = £) (11.20) 


The prediction from this model, and related ones (see Solé 
2003; Solé and Deisboeck 2004; Solé et al. 2008), is that 
unstable tumors close to the critical boundary could be highly 
adaptive but also fragile with respect to small variations of 
instability levels. Although our model is based on fixed parameter 
values, more detailed models involving changing instability levels 
reveal that such levels can increase under Darwinian dynamics, 
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Figure 11.5. Phases in unstable cancer dynamics. The gray area 
indicates the domain where unstable cancer overcompetes the host 
tissue. When close to the separation boundary, small changes of 
instability levels can generate a shift in the behavior (provided 
that normal cells are still present). Two examples of the associated 
potentials are shown, with u = 0.4 and œ = 2.5 (left) and œ = 3.5 
(right), respectively. 


pushing the cancer population close to the critical domain. If 
further instability were introduced (as might be the case while 
using radiation or appropriately designed drugs) the tumor could 
eventually cross the line and decay. 


11.4 Cancer as an Evolving System 


As pointed out by Merlo et al. (2006), cancer is a disease of clonal 
evolution within the body. Although such an idea was formulated 
early (Cairns 1975), modeling and clinical and molecular data 
reveal that evolution and genomic instability in particular, are 
key to carcinogenesis. Theoretical approaches to the origins and 
development of cancer need to take such considerations into 
account. 

The previous models have several important limitations. One 
is their lack of explicit genetic variability and heterogeneity. We 
have already mentioned that tumors are genetically unstable in 
terms of the diversity of cellular populations. This has been shown 
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in many different studies, tracking tumors both in space and 
time (Gonzalez-Garcia et al. 2002; Maley et al. 2006). Since 
spatial patterning is a source of opportunities for diversification 
(Solé and Bascompte 2007), appropriate models should include 
spatial degrees of freedom (by reducing the impact of competitive 
interactions and favoring cooperation) in their ecological de- 
scription of tumors. Moreover, there is also an evolutionary 
element to consider: tumor progression takes place to a large 
extent by exploiting the evolutionary legacy affecting our genomic 
organization. Appropriate models need to consider this feature 
(Spencer et al. 2006). Finally, although our picture of tumor 
dynamics is mainly sustained by the consideration of competitive 
interactions, cooperation might also play an important role in 
driving cancer progression (Axelrod et al. 2007). 
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ECOLOGICAL SHIFTS 


12.1 Change in Ecology 


Ecosystems are complex and adaptive (Milne 1998; Levin 1998; 
Solé and Levin 2002; Pascual 2005), displaying nontrivial forms 
of organization in both space and time. An important component 
of change in ecological systems involves the rapid response of 
some communities to slow changes in external variables (May 
1977; Lefever and Lejaune 1997; Scheffer et al. 2001; Scheffer 
and Carpenter 2003; Rietkert et al. 2004). It has been shown 
in different contexts that constant changes in water availability, 
decline of some plant or fish species or increasing temperatures 
can trigger sudden changes in ecosystem organization. The out- 
come of such changes is often a shift from one stable state to 
another. This can affect a large geographic area, as happened with 
the transition from green to desert in the Sahara, and discussed in 
chapter 3. Such changes are the result of nonlinearities coupling 
external inputs and internal responses. In terms of dynamical sys- 
tems, populations can shift suddenly from one state to another as 
a limiting input parameter changes beyond some given threshold. 

A perfect illustration of the type of sudden changes that 
can occur in complex ecosystems is provided by the dynamics 
of semiarid vegetation. Semiarid and arid habitats represent 
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30 percent of the current ecosystems on our planet, and are 
strongly limited by water availability. Desertification is an increas- 
ing threat to many habitats in the world, having a great impact on 
biodiversity and productivity. Such a process is being accelerated 
by raising human intervention through fires, grazing, habitat loss 
and fragmentation, and decreased water availability. Both studies 
of ancient paleoecosystems (see Foley et al. 2003 and references 
cited) and current arid ecosystems (Kefi et al. 2007; Scanlon 
et al. 2007; Solé 2007) reveal the presence of complex dynamics 
coupling water and vegetation. 

Many other examples of qualitative changes have been 
reported in the literature. Among them, shifts between alternate 
states in shallow lakes provide another compelling example 
(Scheffer et al. 1993). In these systems, the unperturbed, pristine 
state with clear water and submerged vegetation becomes 
replaced by a new state characterized by turbid water (mainly due 
to phytoplankton accumulation). This occurs as a consequence 
of a human-induced increase in eutrophication associated with 
excess nutrient concentration. It was early observed that water 
remains crystal clear under increased nutrient uptake until a 
threshold is passed. Afterward, a shift from clear to turbid occurs 
in an abrupt fashion. Reduction of nutrient availability seldom 
reverses this state, thus suggesting that memory effects (see 
chapter 3) might be at work. 

Another example of ecosystem change leading to extinction 
under the presence of some thresholds is provided by habitat 
loss and fragmentation. Current rates of habitat destruction are 
extremely high. Habitat loss strongly enhances the detrimen- 
tal effects triggered by species introductions, pollution, climate 
change, and hunting. The physical changes associated with habi- 
tat loss and fragmentation include reduction of total area and 
productivity of native areas, isolation of forest remnants, and 
changes in the physical conditions of the remnant fragments. 
These antropogenic changes trigger further community responses 
that sometimes end in a biotic collapse (Wilcove 1987). Multiple 
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imbalances arise, leading to nonlinearities and cascade effects 
throughout the ecological webs. The loss of key species can pro- 
mote subsequent loss of their predators, parasites, and mutualists 
(Terborgh et al. 2001). Given the magnitude and consequences 
of habitat destruction, it is imperative to get enough insight to 
understand the effects of habitat loss on species survival, and 
predict its further consequences. Since economic trade-offs are at 
play, scientists are often faced with the question of how much 
habitat can be destroyed before a certain species goes extinct. 
In order to fully understand the consequences of habitat loss, 
models play a relevant role, particularly in forecasting the effects 
of landscape degradation on web structure and stability. 

Sudden transitions in ecosystems are likely to be common. 
Such tipping points represent an important problem, not only 
in the context of ecology but also in other fields, including 
medicine, geology, and finance (Scheffer 2009). One obvious 
question is, how predictable are they: can we detect warning 
signals? (Scheffer et al. 2009). Below we consider two examples 
of critical transitions in arid ecosystems and outline the problem 
of early signal detection. 


12.2 Green-Desert Transitions 


We choose here a very simple type of model that has been intro- 
duced in two different forms in the literature (Klausmeier 1999; 
Shnerb et al. 2003). They are both examples of the nonlinear 
interactions between vegetation and water availability and they 
were both introduced as spatially extended interactions. These 
systems are known to exhibit a variety of spatial arrangements of 
vegetation resulting from a self-organization process that has been 
successfully modeled.’ Some examples of such vegetation patterns 


! Specifically, these systems seem to exhibit a so-called Turing-like instability, 
where nonlinearities together with spatial processes such as diffusion can 
generate the observed regular structures. 
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Figure 12.1. Spatial patterns of vegetation cover are known to 
spontaneously develop under a wide range of conditions in ecosys- 
tems. In (a-c) we show three examples of such patterns in (a) 
arid, (b) semiarid, and (c) Mediterranean ecosystems (adapted rom 
from Shnerb et al. 2003). These three snapshots involve increasing 
rainfall levels. The basic interactions between water availability and 
vegetation cover are summarized in (d). 


are shown in figure 12.1a—c for arid, semiarid, and Mediterranean 
ecosystems, respectively. 

Without taking into account this spatial component, we can 
consider the mean field model just by looking at how water 
and vegetation cover interact. The basic scheme of plant-water 
interactions is shown in figure 12.1d. Here w and 7 indicate 
water and vegetation cover, respectively. Water is an external 
input (rain) that we indicate by means of a constant œ, and 
it is depleted by run-off proportional to the current amount 
(—ew). Additionally, we need to introduce the consumption 
of water by vegetation. Such consumption can be linear: the 
higher the plant density, the larger the water consumption. But 
it can also be nonlinear, represented here by a quadratic term. 
Such a term tries to capture underlying plant-plant interac- 
tions, in particular possible facilitation effects: as more plants 
become present, better local conditions enhance further growth. 
Moreover, plants die at some density-independent rate m and 
grow by exploiting available water with some given efficiency 
p. These basic assumptions allow a two-dimensional dynamical 
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system: 
dw 
— =a-—ew—vwI(n) (12.1) 
dt 
ae = —mn + pwl (n) (12.2) 
dt 


where the function I introduces the functional form of the 
water-vegetation interaction. Here the two cases are considered: 
T(n) =n (Shnerb et al. 2003) and the quadratic dependence asso- 
ciated to facilitation effects T(z) = n? (Klausmeier 1999). The 
latter can easily be interpreted: the “encounter” of neighboring 
plants in a given area enhances the likelihood of each to exploit 
available water resources. For example, a given plant with well- 
established roots will help other plants to establish themselves. 
This model can be further reduced by considering that water 
dynamics proceeds much faster than vegetation dynamics. On a 
first approximation we would have dw/dt ~ 0 and thus we can 
replace w in the second equation by 
a 


we ITO) (12.3) 


which gives, once introduced in dn/dt, a one-dimensional 
system: 

dn apr (n) 

@ sio 


This general equation contains two different cases, associated 


(12.4) 


to a linear and nonlinear T (7) term. It also defines two different 
types of transitions and will be analysed below. 


12.3 Continuous Transitions 


Let us first consider the linear case. Here the dynamics is de- 
scribed by: 
dn aon 


ee ae = 12. 
dt €+n il Me) 
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Figure 12.2. Transitions from vegetated to desert states in Shnerb’s 
et al. model. In (a) we plot the bifurcation diagram for € = 1. The 
corresponding parameter space is displayed in (b) where the two basic 
domains are indicated, separated by the critical boundary jz, = €. 


and it is easy to show that two fixed points are possible and 
correspond to the desert state n* = 0 and the vegetated state 
nt = e (12.6) 
m 
The stability of these two points can be determined from A(0) = 
pa/é — m, which gives us a critical condition for the transition 
from one phase to the other. 

Using the notation = pa/m, we have a desert state for 
jL>€ =p, and a vegetated state otherwise. This jz (control) 
parameter is directly dependent on water availability and plant 
efficiency, whereas it is inversely proportional to vegetation sur- 
vival rate. It thus properly incoporates all the relevant trade-offs 
introduced by the model. In figure 12.2a we show the bifurcation 
diagram for this system (when €=1) and the corresponding 
phase space, where the critical line u, = € separates both scenarios 
(figure 12.2b). As we can see, the transition is second-order, with 
a continuous change of the order parameter 7”. 
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12.4 Sudden Shifts 


2 is rather 


The behavior for the second case, when T(n) = n 
different. In this case, three fixed points are at work. The first is 
again the trivial one, n* = 0. This fixed point is always stable, 
since A,,(0) = —1 < 0, indicating that small vegetation covers 
will end up in a desert state. The two other fixed points are 


obtained from the quadratic equation: 
n? —un*+e=0 (12.7) 


and this gives 


t= su + y u? — áe] (12.8) 


which exist provided that 
U > Me = Dye: (12.9) 


The nontrivial fixed points (when they exist) are separated 
by a gap from the zero state. Since the latter is always stable, 
for consistency the two other fixed points, n*, and nł, will be 
unstable and stable, respectively.” The corresponding bifurcation 
diagram and phase plot appears in figure 12.3a. This diagram 
reveals something of great importance related to the potential 
irreversibility of the shift. Starting from a vegetated state, as 
we change u (say, water availability) below u., the transition 
toward desertification cannot be reversed by simply moving back 
to u > He. The system is now locked in the new attractor and 
no jump back will occur, unless we enter again in to the domain 
n > në. 

We can show the emergence of a bifurcation leading to a three- 
fixed-point situation by plotting the two components defining 


2 Since n* > n* > n*, they define an ordered sequence on the real line. The 
all is Se eae q 
stability of n% implies that trajectories come out from 7* and similarly flow into 
ais which is thus stable. 
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Figure 12.3. Transitions from vegetated (green) to desert states in 
Klausmeier’s model. In (a) we plot the bifurcation diagram for € = 1. 
The corresponding parameter space is displayed in (b) where the 
two basic domains are indicated, separated by the critical boundary 


Ue =e. 


our reaction terms. To make the analysis simpler, let us fix € = 1 

and divide all terms in the differential equation by m. In this way 

we have a single parameter model: 
dn pun 
dt 1+% 


where time has been rescaled by the mortality rate m. Since 


—n (12.10) 


n* = 0 is already a fixed point, the two functions are reduced 
to a nonlinear one, namely 

un 
1+? 


and a constant one, y2(a)= 1. The nontrivial fixed points are 


yi(n) = (12.11) 


reached at the intersection y; (7) = y2(2), which means the cross- 
ings between the horizontal line y = 1 and the curve yp. 

A simple analysis reveals that the critical value u, = 2 is indeed 
the bifurcation point of our system (figure 12.4a). The curve y2 
is such that y2(0) = 0, but y2 —> 0 for n — oo. Since it is 
positive for the n > 0 meaningful interval, we should expect to 
find a maximum at some given value 7,,. The extrema of y2() are 


Vegetation cover n 


Figure 12.4. The presence and character of the fixed points in 
Klausmeier’s model can be determined by using a simple graphical 
method (a) where the two components of the dynamics are plotted 
(see text). In (b) we show the predicted potential function ®, (7) for 
different values of the bifurcation parameter jz. As we cross p, = 1, 
the equilibrium represented by the marble shifts from a nonzero 
(green) vegetation to the desert state. 


obtained from the condition dy,/dn = 0, which gives n,, = 1. 
Now we can calculate the value of the maximum y,, taken by our 
function at this point. We obtain 


(2m) = 5 (12.12) 


and this immediately allows us to determine the conditions for 
crossing with yı = 1: if u < 2 no crossing is possible since the 
curve is below y = 1 (figure 12.4a). Once we cross this value, two 
intersections are expected to happen. 

Finally, we can also calculate the form of the potential function 
®,,(z) whose minima are associated with desert or vegetated 
states. By using Klausmeier’s model, we obtain: 


2 
,(n) = > — p (n — arctan n) (12.13) 


which is depicted in figure 12.4b for different values of jz. As 
predicted by our previous stability analysis, the marble at the 
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bottom of ®,,(7) will move from n* > 0 to the desert state n* = 0 
as we cross the critical level ue. 


12.5 Metapopulation Dynamics 


Let us now consider the problem of species extinctions under 
habitat loss. The simplest, first approximation to this problem 
is obtained by exploring the consequences of habitat loss in a 
metapopulation (Hanski 1999; Bascompte and Solé 1996). A 
metapopulation can be defined as a set of geographically distinct 
local populations maintained by a dynamical balance between 
colonization and extinction events. Such events can be described 
in terms of state transitions happening with a given probability. 
Let us start this section by revisiting the Levin’s (1969) model, 
which captures the global dynamics of a metapopulation. 


The basic model reads: 


dx 
dt 


where x is the fraction of patches occupied, and c and e are the 


=cx(1—.x) -—ex (12.14) 


colonization and extinction rates, respectively. This model has an 
extinction state xj = 0 and a nontrivial solution given by Levin’s 
steady population, x7 = 1 — e/c. The colonization rate has to be 
larger than the extinction rate for the metapopulation to persist. 
The resulting model has the same structure as the one already seen 
in chapter 9 on epidemics. It can also be rewritten as a logistic-like 
model as follows: 
as =(c—e)x (1 — =) (12.15) 
dt x$ 
where the condition for stability of x7 becomes obvious since we 
need c > e in order to have a positive growth rate. 
One can easily introduce habitat loss into the framework of 
model (12.14). If a fraction D of sites is permanently destroyed, 
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Figure 12.5. Phase changes in metapopulations under habitat loss. 
In (a) the two phases of Levin’s model, namely extinction (white 
area) and stable (gray) are shown. In (b) the corresponding order 
parameter x*(c, D) is shown (for € = 0.1.) Either if colonization is 
low or habitat destruction too high, extinction occurs by crossing the 
critical boundary D, = 1 — e/c. 


this reduces the fraction of vacant sites that can potentially be 
occupied. Model (12.14) becomes now: 
dx 


— =cx(l1— D — x) — ex. (12.16) 
dt 


This model has two equilibrium points (to be obtained from 
dx/dt = 0): x* = 0 (extinct population) and x* = 1— D—e/c. 
The latter decays linearly with habitat loss, becoming zero when 
1 — D — e/c = 0. This condition gives a critical destruction 
level D, = 1 — e/c indicating that a nontrivial dependence exists 
between available patches and the species-specific extinction and 
colonization properties. Once we cross this threshold (i.e., if 
D > D,) the population goes inevitably extinct. This situation is 
illustrated in figure 12.5. Here the critical line separates the two 
qualitative types of behavior. As we approach the critical value D, 
by increasing the amount of habitat destroyed, the frequency of 
populated patches decays linearly, becoming zero at the boundary. 
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This process is actually accelerated by considering local dispersal 
(Bascompte and Solé 1996). 

The previous models can be generalized in order to observe 
first-order transitions by incorporating an additional constraint. 
It explains a common observation from the analysis of field data: 
it appears that species occur, at any one time, either in most 
potential sites or only in a few. This bimodal pattern can be 
easily explained by means of the so-called Allee effect, which 
introduces a negative contribution to population growth at low 
densities (Stephens and Sutherland 1999; Fowler and Ruston 
2002; Liebhod and Bascompte 2003; Courchamp et al. 2009). 
This can be for a number of reasons, from mate shortage, failure 
to exploit available resources, to lack of (required) cooperation. It 
can also be induced by anthropogenetic means (Courchamp et al. 
2006). 

The Allee effect can be easily incorporated within Levin’s 
formulation by assuming that there is some additional constraint 
limiting the survival of local populations. This can be done by 
properly modifying the model (Dennis 1989; Amarasekare 1998; 
Keitt et al. 2001; Courchamp et al. 2009): 

a a (:-=) (12.17) 


j. * * 
dt x7 xF 


Here a new parameter 4, such that 


e 
0<a<1-- (12.18) 
c 


and defining a threshold fraction of occupied sites below which 
populations go extinct. The new model involves a cubic equation, 
and thus three fixed points are present. Two of them are identical 
to those found in the original model, whereas the third, x* = a, 
is new. It introduces an intermediate state that strongly changes 
the structure of the bifurcation diagram. 

It is easy to show now that the extinction point x* = 0 is always 
stable, in sharp contrast with the original Levin’s model, where 
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Figure 12.6. Levin’s model with multiple states due to the Allee 
effect. Here we use a = 0.3 and c = 0.7 and plot the bifurcation 
diagram (a) and the paramater space (b) for different values of 
extinction. As we can see in (a) a bistable state is generated, the 
x* = a fixed point the unstable one. 


this state depends on the parameters. By studying the other two 
fixed points we can see that the x7 is now stable if c(1 — a) < e 
whereas the second is unstable (it does not exist outside this 
domain). The resulting bifurcation diagram is shown in figure 
12.6a. We can see that it describes a first-order transition: once 
we cross the e, = c(1 — a) critical boundary, the population 
shifts to zero. The model can be further explored (Amarasekare 
1998) by incorporating the destruction rate: 

Z = (c(l — D — x) — e)x 5 = 7 


We leave to the reader the analysis of this extended model. Let us 


) (12.19) 


just mention that metapopulation persistence now requires more 
strict conditions. 


12.6 Warning Signals 


If critical changes can generate catastrophic events by slowly 
tuning a key parameter around the transition point, it would be 
desirable to have warning signals indicating that such changes are 
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likely to occur (Scheffer 2009; Scheffer et al. 2009). Some of such 
early warning signals have been shown to exist. They are associ- 
ated to both time- and spatial-dependent ingredients. Specifically, 
it has been shown that critical slowing down (see section 3.4) 
can be expected when approaching the threshold (Chrisholm and 
Filotas 2009). Here natural fluctuations would be amplified by 
the internal nonlinear dynamics and a high variance of population 
change should observable. Similarly, some special types of spatial 
patterns are known to arise close to the transition (Riekert et al. 
2004; Kefi et al. 2007; Kefi 2009). As it occurs with magnetic sys- 
tems (chapter 1), spatial correlations experience dramatic changes 
close to criticality. It has been reported from both ecological and 
climate studies that marked trends in variance measures typically 
occur close to criticality (Lenton 2008; Anderson 2009; Dakos 
et al. 2008; van Nes and Scheffer 2007). 

A special problem, addressed in terms of phase transitions 
in ecology, concerns the propagation of fires. Forest fires are 
known to be an active component of ecosystem dynamics and 
are an important evolutionary force (Bond and Keeley 2005). 
They promote and maintain diversity by creating a range of 
opportunities in local habitats of varying conditions. Fires are 
actually similar to herbivores (Bond and Keeley 2005), and, 
as mentioned in chapter 9, their dynamical behavior is closely 
related to that exhibited by epidemic spreading. Although fires 
are a natural outcome of ecosystem dynamics, their frequency 
and size are growing across the tropics. Their interactions with 
increased habitat degradation, logging, and climate change have 
been responsible for huge carbon emissions (Cochrane 2003). By 
analyzing the patterns displayed by Amazonian forest fires, it is 
possible to predict some potential scenarios of megafire spreading 
(Pueyo 2007; Pueyo et al. 2010). These studies suggest that, if 
the most adverse forecasts are realized, the second-order transition 
expected from an epidemics-like behavior would be replaced by 
sudden (first-order) changes. 


13 


INFORMATION AND TRAFFIC JAMS 


13.1 Internet and Computer Networks 


One seemingly general property of transport networks, including 
roads, Internet, the power grid, or a highway, is the undesired 
emergence of jamming and congestion. As the flow of matter, 
energy, and information increase beyond a given threshold, 
problems start to pile up. This situation is common in roads 
when many vehicles become jammed, but it is also a problem for 
the power grid when energy demands exceed the grid’s capacity. 
Congestion is actually one of the oldest examples of a phase 
transition that we can think about. Who hasn’t experienced a 
traffic jam or an Internet storm? But few people realize that, both 
in highways and computer networks, it involves a phase transition 
between fluid behavior and congestion (Nagel and Schreckenger 
1992; Nagel and Rasmussen 1994; Huberman and Lukose 1997; 
Ohira and Sawatari 1998; Solé and Valverde 2000). The presence 
of two phases has been reported in different field studies. An 
example is shown in figure 13.1, where the fluctuations of packets 
flowing through a given link in a computer network are displayed. 
Over some intervals, the link becomes congested by high traffic, 
which saturates it, whereas for other intervals a lower density 
fluctuation is seem (Takayasu et al. 2000). 
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Figure 13.1. An example of time sequence of flow density fluctu- 
ations observed at intervals of 0.1 s in an Internet node (redrawn 
from Takayasu et al. 2000). An intermediate interval reveals a very 
high flow of traffic, making the system highly congested. The other 
two intervals reveal a highly fluctuating, but otherwise fluid, traffic 
regime. 


Parallel multiprocessor networks have been shown to display 
complex dynamics and a phase transition separating a congested 
from a noncongested phase, both in the real and the simulated 
(Ohira and Sawatari 1998; Solé and Valverde 2000) counterparts. 

In real systems, there is a coupling between the system’s 
dynamics (exhibiting fluctuations) and the human agents, which 
both create and respond to such changes. This feedback is 
clear if we take into account how users interact within the 
Internet while searching for information. Each user sends and 
receives information and as a consequence contributes to the 
traffic. However, when the flow becomes too slow because of a 
congestion event, the behavior of the users changes: congestion 
reduces their pressure on the system and eventually a more fluid 
state is recovered. Such feedback has been introduced in models 
of Internet traffic (Valverde and Solé 2002, 2004), which have 
been shown to consistently reproduce several key properties of real 
computer traffic, suggesting that the complex fluctuations might 
result from a critical state. 
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In its simplest form (to be considered here, see figure 13.2a), 
we can imagine a computer network as a set of connected 
elements, some of which (the so-called hosts) create and send 
packets to some other chosen host node. The other type of 
element, the so-called routers only canalize the traffic through the 
web and cannot store packets. Instead, hosts can store arriving 
packets in a queue, which is constantly shortened as the host 
processes the incoming messages (which are sequentially piled 
following their order of arrival). Once the packet is received, it 
is either removed (if the queue is empty) or appended at the 
tail of the host queue. These hosts are assumed to be randomly 
scattered over the lattice. Such a pattern matches some special 
types of well-known computer architectures (Hillis 1984) and has 
been explicitely implemented in a grid structure (Bolding et al. 
1997). At low rates of packet injection, we should expect fluid 
traffic, whereas for high rates, the system will inevitably become 
congested. Below we develop the simplest mean field approach 
for this problem, showing that two phases separated by a sharp 
transition should occur. 


13.2 Mean Field Model of Traffic Flow 


Although the grid shown in figure 13.2a implies spatial cor- 
relations between nodes, we ignore them and assume that the 
nodes are somewhat mixed, displaying an average number of 
connections (k) among them. We also consider queues having 
a potentially unlimited length. A simple mean field model can be 
obtained for the total density of packets, namely T (£) = N(t)/L?. 

The number of traveling packets increases as a consequence 
of the constant pumping from the hosts, which are present at a 
constant density p. If each host sends packets at a rate u the total 
rate of packet emission is pA. However, packets are efficiently 
removed from the system if the lattice is not too congested (i.e., 
if free space for movement is available) but accumulate as a 
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Figure 13.2. A simplified model of computer traffic. Here (a) the 
elements are arranged on a square grid with L x Z nodes. Two 
types of nodes are considered: hosts (indicated with filled circles) 
and routers (plain squares). Hosts send and receive packets, whereas 
routers simply redirect them through the grid. The grid is folded in 
such a way that we form a torus (b). In a model of computer traffic, 
different nodes receive different amounts of traffic load, and different 
levels of congestion will be observed. In (b) this is indicated by lighter 
gray values. 


consequence of already jammed nodes. We can roughly state that 
once the number of packets exceeds the number of lattice sites, 
congestion will lead to packet accumulation. Since the total lattice 
size is L?, we can consider this value the critical limit of packet 
load that the system can manage without getting congested. 
The time evolution of the density of packets will follow the 
mean field equation for the number of packets N: 
= = pul? -atn (1 — *) (13.1) 
As we can see here, we adopt a basic criterion for congestion: if 
the number of packets in the system is smaller than the grid size 
(N < L?), no congestion takes place, whereas if N > L?, under 
our mean field picture, it is likely that congested traffic will be 


present. 
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Using the previous definitions the model can be normalized 
and it reads: 
dT 


— = pu -a(r (1 -T) (13.2) 
dt 


The last term on the right-hand side indicates the rate of removal, 
which is proportional to the number of input pathways (number 
of neighbors (&)) available to incoming packets and the efficiency 
a with which they are removed. The rate œ corresponds to the 
inverse of the average lifetime of packets (Fuks and Lawniczak 


1999). 


13.3 Fluid Flow and Congestion Explosion 


In our approximation, we consider a fixed, constant injection of 
packets from the hosts. For A constant (i.e., when no reactive 
responses are assumed to exist) we have the two fixed points, 


namely: 
1 4p 
Te =- }|l+,/1-—~ 13. 
thal] as 
These two fixed points exist provided that u < u, where 
k 
p=) (13.4) 
4p 


The stability of these two fixed points can be determined from 
the sign of A,, = —a (k) [1 — 2r*], which in our case reads 


dy (T's) = a (k) h +2,/1- z| (13.5) 


he 


and thus we conclude that F+} is unstable whereas the lower 
branch T— is stable. Outside this parameter domain, the system 
undergoes an explosive growth of congestion. 

The corresponding bifurcation diagram is shown in figure 
13.3a. This is a rather special diagram, involving a single 
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Figure 13.3. (a) Bifurcation diagram for the traffic congestion 
model. Here the bifurcation parameter is the rate of packet emission 
(up) by each element in the grid (we fixed (k) = 2 and p = 0.1). 
For u < u, two fixed points are present but only one of them (the 
lower branch) is stable. After the critical emission rate is crossed, the 
system diverges and congestion increases indefinitely, leading to the 
system’s collapse. In (b) an alternate representation is given, showing 
the two basins of attraction defined by the system dynamics. The gray 
zone involves all the initial values that lead to the stable, fluid state. 


equilibrium state toward which a range of initial conditions 
converges. The two branches obtained from our analysis are 
indicated for y < u, and we can see that the diagram is similar 
to other plots associated with first-order transitions. But here 
the lower branch is not the unstable one: the order has been 
exchanged, meaning that initial conditions starting from values 
r(0) > r4 will lead to increasing congestion. Similarly, the lack 
of a (finite) fixed point for u > u, needs to be interpreted. But 
the solution is simple: once we have enough packets in our system, 
a runaway effect generates an explosive increase that creates a 
density of packets crossing the F = 1 value. Afterward, the last 
term in the mean field differential equation will turn positive and 
congestion will allways rise. The two-phase behavior is illustrated 
by the trajectories following the system, as shown in figure 13.4a. 
The corresponding phases are also shown in figure 13.4b, sepa- 
rated by the critical line u, = a(k)/4p. 
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Figure 13.4. Phase transition in the mean field model of traffic. In 
(a) we show the trajectories obtained numerically using as initial 
condition r (0) = 0.5 and different values of u (other parameters 
are (k) = 2 and p = 0.1), which gives a critical rate u, = 5. For 
u values below critical, the system’s activity relaxes to a stable value 
T* (see text). In (b) we show the phase space for this model using 
u and p as control parameters. The shadowed area corresponds to a 
divergent number of packets. 


13.4 Self-Regulated Traffic 


In standard computer networks, congestion levels are bounded 
for two good reasons. First, queues are not infinite and thus 
unbounded growth of stored packets cannot occur. Second, 
and most important, the rate of packet emissions cannot be 
independent of congestion levels. The feedback between the two 
components can trigger complex fluctuations (not taken into 
account here) and scaling behavior (Huberman and Lukose 1997; 
Valverde and Solé 2002, 2004). In this context, it has been 
shown that model computer networks spontaneously evolve to 
an intermediate state close to the critical boundary separating 
fluid traffic from congestion. Such a situation is known as self- 
organized critical and has been suggested to be present in many 
different situations both in nature and society (Bak 1997). 


13.4 Self-Regulated Traffic 155 


Under our mean field approximation, a simple (but important) 
modification allows us to show that such a critical state is 
achieved. Considering that the level of packet occupation 1 — F 
affects not only packet dynamics but also rate of emission, we can 
indicate a tentative model as follows: 

dT 


a (ou — a(k)T) (1 —T) (13.6) 
t 


which now defines two stable equilibrium points, namely T = 1 
and T = pyu/a. We can easily check that 


ar ar 
(=) CS] =r . A o laa mete 


and thus a nonsaturated point will be stable provided that 
U < u: =a(k)/p. Otherwise, congestion rises until complete 
collapse occurs. This result is consistent with a continuous phase 
transition, and the critical point is clearly similar (up to some 
constant factor) to the one just discussed above. Now however, for 
a fixed u < Ue, a single solution exists, where the system achieves 
a steady state with intermediate levels of traffic. 

As we have mentioned in other situations, the architecture 
of the actual web of computer interactions will also have an 
impact on its behavior. The Internet is known to be a complex 
heterogeneous network, with most nodes having only a few 
connections whereas a handful of nodes display a huge number of 
links (Yook et al. 2002). Such a topological pattern makes these 
webs very efficient but fragile, and the strategies of packet routing 
play a relevant role in making traffic efficient. In this context, 
other types of phase transitions can be found. An interesting 
example concerns the depth of routing tables (Valverde and Solé 
2004). As we outlined above in our lattice model, routers redirect 
packets toward their destination. However, in a multimillion-site 
web the routing tables are limited to a maximum depth. Beyond 
this horizon, the packet would travel randomly. Although it 
might appear optimal to have huge routing tables spanning the 
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whole network, it can actually be shown that there is a critical 
path horizon beyond which packets are successfully delivered.' 
This critical horizon corresponds to a phase transition (Valverde 
and Solé 2004) and shows that an appropriate combination of 
routing and randomness leads to an optimal information flow at 
a low cost. 


! In fact, a system incorporating full tables will collapse, since some particular 
links will be typically chosen as key elements within the shortest paths and 
experience traffic jams. 


14 


COLLECTIVE INTELLIGENCE 


14.1 Swarm Behavior 


One of the most remarkable events in the evolution of life was the 
emergence of social behavior among insects. The generation of 
complex societies of cooperating individuals represented a major 
leap in the history of our biosphere. It is difficult not to become 
fascinated by the plethora of patterns displayed by the collective 
work of ants, termites, bees, and social wasps. The huge nests 
created by termites or the army ant raid patterns traveling through 
the rainforest are just two examples (Holldobler and Wilson 
2009). The ecological relevance of these insects is illustrated by 
the following fact: in some rainforests, the dry weight of ants 
and termites is about four times that of all the other land animals 
(Holldobler and Wilson 1990). 

Social insects build complex structures, scan their environ- 
ments in search of available resources, perform group-based 
decisions (Deneubourg and Goss 1989; Camazine et al. 2002), 
and—in many ways—remind us how brains work (Hofstadter 
1980; Solé and Goodwin 2001). In social insects, while colonies 
behave in complex ways, the capacities of individuals are relatively 
limited. But then, how do social insects reach such remarkable 
ends? The answer comes to a large extent from self-organization: 
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Figure 14.1. Symmetry-breaking in ant colonies. Here two exam- 
ples of experiments, illustrating how the collective behavior of the 
colony leads to a broken symmetry. In (a) escaping ants choose 
preferentially one of the two possible exits (symmetrically located on 
both sides of a circular arena). The image in (b) shows ants moving 
through one of two bridges between the nest and a food source. 
After a while, all use one path (picture kindly provided by Ernesto 
Althsuler [a] and by Guy Theraulaz[b]). 


insect societies share basic dynamic properties with other 
complex systems (Millonas 1993; Gordon 1999, 2010; Detrain 
and Deneubourg 2006; Sumpter 2006) and exploit bifurcations 
and nonlinearities as their source of internal organization and 
decision making (Deneubourg et al. 1989; Bonabeau 1996). 
Moreover, some trends displayed by the global organization of 
insect societies, such as the diversity of morphological casts, nest 
organization, and hierarchy have been related to phase transition 
phenomena (Oster and Wilson 1978; Valverde et al. 2009). 
Several experiments clearly indicate that amplification of per- 
turbations by means of individual interactions leads to collective 
responses. Two examples are shown in figure 14.1. In the first (a) a 
small group of ants in a circular arena suddenly experience a panic 
reaction to an externally introduced chemical signal. The arena 
is enclosed by a plastic tube with two identical exists, but ants 
tend to choose one of them: a symmetry-breaking phenomenon 
pervades the panic behavioral response. The second example, to 
be analyzed here (b) deals with how ant colonies can choose 
among two given bridges connecting their nest with a food source. 
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Foragers follow the chemical trails left by their fellow ants. The 
more pheromone is deposited on a given path, the higher the 
probability that ants follow that path and leave further chemical 
marks. When the paths have different lengths, ants typically end 
up choosing the shortest, thus solving an optimization problem 
(Deneubourg et al. 1983, 1987; Millonas 1992; Bonabeau 1996; 
Beckers et al. 1993; Deneubourg et al. 1990). But even when the 
two bridges are identical, the amplification process can lead to a 
choice, thus breaking the symmetry of the system. Here we will 
consider this special case. 

In this chapter we study two different types of phase transi- 
tion phenomena that have been reported from the experimental 
analysis of social insects. They illustrate how collective behavior 
can allow colonies to make decisions and display self-organization 
behavior. 


14.2 Foraging and Broken Symmetry 


The two-bridges problem can be formalized as follows. Let us 
indicate by xı and x2 the concentrations of trail pheromone in 
each branch. These concentrations will be somewhat proportional 
to the number of ants exploring each path. The equations for 
these quantities will read (Nicolis and Denebourg 1999): 

a = uqi P(x, x2) — Vx] a = uq P(x, x2) — Vx 


dt 
(14.1) 


where u is the rate of ants entering each branch from the nest, 
qi the rate of pheromone deposition in the 7 — th branch, and 
v the rate of pheromone evaporation. The functions P;(x1, x2) 
introduce the choices made by ants (which bridge is chosen) given 
the pheromone concentration. A possible response function is 
described by (Beckers et al. 1992; Deneubourg et al. 1990): 
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Figure 14.2. Symmetry breaking in the food search experiment. 
When two identical bridges are used, the amplification process 
associated to pheromone fields forces the collective to make a choice, 
breaking the initial equivalence. 


where O(x1, x2) = 2 ja1,2(%; + K)* andi = 1,2. It can be 
shown that this model allows interpretation of the minimization 
process in terms of a symmetry-breaking bifurcation. 

Let us consider first the symmetric case where q) = q2 = q 
and thus the set of equations: 


dxı — (xi + Ky dx — (x2+ Ky’ E 


dt E ma O(x1, x2) 7 l “dt ta 


(14.3) 


dt ~ T Ola 


This symmetric scenario is not very interesting in terms of 
decision making, since there is no better solution. But the exercice 
is useful and shows that a choice can be made by amplification of 
fluctuations. The same principle pervades the more general case, 
to be explored later. It is not difficult to show that one possible 
solution to this system is the symmetric fixed point xf = x} = x* 
when ants equally distribute themselves over both branches. For 
this special case, we have P;(x1, x2) = P(x) = uq/2 and thus 
a single equation dx/dt = uq/2 — vx that gives a fixed point 
x* = ytg/(2v). The second state would correspond to the choice 
of one of the branches. Since the total pheromone concentration 


2.0 
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is 2x* = uq /v, we can see that 


2 
(2 -x)@+Kpex(-x+K) (14.4) 
v v 

After some algebra, this gives the new fixed points x = (xf, x3) 

with 
»_ #4 ugy? x Hq ugy? 
= it Pe), os ee a (EL) a Ke 
“a 2v H a, *2 2v Ge, 
(14.5) 


and x; = uq /v — x}, (i # j). This pair of fixed points will exist 
provided that ug /2v > K, which allows us to obtain the critical 


point at 


2K 
ka (14.6) 
q 


Below this value, the only fixed point is the symmetric case with 
identical flows of ants in each branch. This critical curve allows 
to define two phases, as indicated in figure 14.2c, where we 
plot u: = u. (q) using K = v = 1. The gray area marks the phase 
where collective choices are possible, whereas the white area 
is associated with random behavior and no amplification of 
fluctuations. In this latter domain, no symmetry-breaking is 
possible. Since the symmetric point x* becomes unstable and 
we have two new symmetrically located branches around it, 
these two new solutions must be stable, provided that u > ue. 
The corresponding bifurcation diagram is shown in figure 14.3. 
This symmetric model can be generalized to (more interesting) 
asymmetric scenarios where the two potential choices are different 
(see Detrain and Deneubourg 2009 and references cited) either 
because the food sources have different size or because paths have 
different lengths and the shortest path must be chosen. 

Transitions in colony behavior associated to food search have 
been characterized in other contexts, sometimes displaying first- 
order transitions (Beckman et al. 2001; Loengarov and Tedeshko 
2008). 
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Figure 14.3. Bifurcations in the symmetric bridge model of collec- 
tive choice. Here we have used K = 1, q = 2, v = 1 and tuned the 
parameter p (here p, = 1). The straight line gives the symmetric 
state x* = uq /2v. The two branches define the two alternative states 
where one branch has been chosen. The corresponding patterns of 
ants walking on each branch are indicated at right. 


14.3 Order-Disorder Transitions in Ant Colonies 


In this section we will address a different type of dynamical 
pattern exhibited by ant colonies, namely the presence of global 
activity oscillations (Franks and Bryant 1987; Cole 1991a). This 
type of phenomenon has been described in different genera of 
ants, but has been particularly well studied using colonies of the 
genus Leptothorax. Briefly, these small colonies (composed of 
just fifty to one hundred individuals) show a remarkable time- 
dependent pattern of activity, where periods of no movement are 
followed by bursts of all-active individuals. These bursts are the 
outcome of local interactions among ants and appear at almost 
periodic intervals over time (figure 14.4). 

By using groups of ants of different sizes, a series of experi- 
ments revealed that the density p of ants enclosed within a given 
arena had a great impact on the observed dynamics. Specifically, 
it was found that the pattern of synchronization exhibited by the 
colony was not present at the individual level. In other words, the 
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Figure 14.4. Fluctuations in colony activity are known to occur in 
many species of ants, such as those of the genus Leptothorax. These 
ants (a) form small colonies (picture by Bernhard Seifert) exhibiting 
almost regular bursts of global activity (b) characterized by intervals 
of little activity followed by intervals where most or all ants are 
performing some task (adapted from Cole 1991a). 


global clock was not the result of the synchronization of multiple 
individual clocks, as it occurs in other situations (Strogatz 2003). 
Coherent bursts of activity appear close to a given critical density 
Pe: but there is no regular behavior below the threshold. Instead, 
individual ants appear to display chaotic behavior (Cole 1991b). 
The coherent pattern displayed by the colony as a whole has been 
explained in terms of nonlinear dynamics (Goss and Deneubourg 
1988; Solé et al. 1993; Cole 1993; Boi et al. 1998; Delgado and 
Solé 1997, 2000). In an inactive colony, an individual can become 
spontaneously active, moving around and exciting neighboring 
ants, which also move and can activate others. A percolation-like 
phenomenon takes place at the critical density. 

The simplest mean field model that can be formulated is 
closely related to the one obtained in chapter 9 for the propaga- 
tion of epidemics. As suggested in Cole (1991c), the propagation 
of activity in ant colonies is similar to the epidemic-spreading 
model. Let us now consider our population composed. of two 
groups, namely active and inactive ants, whose numbers will 
be indicated by A and 7, respectively. Using the normalized 
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values x = A/N and y = Z/N, and assuming that the density 
of occupied space is p € [0, 1], the mean field equation for the 
fraction of active ants is: 
aia =ax(p —x)— yx (14.7) 
dt 
An important difference in relation to epidemic spreading is that 
the density of ants is now a key control parameter. 
The equation displays two equilibrium points. The first corre- 
sponds to no activity at all (x = 0) while the second is associated 
to the presence of a given level of activity, namely 


pape (14.8) 
q 
The colony will display nonzero activity (i.e., xï will be stable) 


provided that 
A, (0) =ap—y >0 (14.9) 


which, in terms of the density of ants corresponds to the condition 
p > pe =y/a. 

The deterministic model predicts that no activity will be 
present at subcritical densities, whereas experiments indicate that 
activity is present, largely due to random individual activations, 
at small densities. Moreover, it has been found that fluctuations 
increase close to the critical point, consistent with real data (Solé 
et al. 1993; Miramontes 1995). This is illustrated by the shape of 
the associated potential, which now reads: 


3 2 
p(x) = a — (ap — > (14.10) 


Three examples of this potential are shown in figure 14.5c for 
subcritical, near-critical, and supercritical cases. As expected, the 
potential flattens as criticality is approached. 

The observed variance in colony fluctuations is a consequence 
of critical slowing down (see chapter 1), and for this model it can 
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Figure 14.5. (a) The two phases displayed by the mean field model 
on the (œ, p) paramater space. The gray area corresponds to the 
active colony state, where colony activity x* > 0 is present. In (b) 
the corresponding bifurcation diagram is shown (in both cases we fix 
y = 0.2). Three examples of the associated potential are shown in 
(c) for p = 0.3, 0.55 and p = 0.7 (from top to bottom). 


be shown that the relaxation time t(p), scales as 


x(T(p)) 
T(p) = ~ | ss (= p) (14.11) 
æ Jeo xP = pe — x) 

As we approach criticality, random activations are likely to prop- 
agate for a while before global activation dies out, thus providing 
the source for higher activity levels and broader fluctuations. 

It seems strange that ants display global colony oscillations. 
A simpler pattern of activity, where all ants are allways active at 
some low level, would appear more natural. In both cases, the 
same global average activity could be obtained. There is, however, 
a good reason for the oscillations. A detailed study of the two 
alternatives (Delgado and Solé 2000) shows that, functionally, 
a fluctuating system actually outperforms the alternative (not 
observed) constant-activity scenario. 


14.4 Structure and Evolution in Social Insects 


We have presented examples illustrating two ways of solving 
an optimization problem based on collective behavior. How- 
ever, insect colonies display multiple levels of complexity, and 
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the presence of critical thresholds is widespread. Thresholds in 
behavioral responses have been identified in the transition to 
social behavior, the distribution of workers into castes, the emer- 
gence of hierarchical organization in wasps, and aggregation phe- 
nomena. Moreover, some key traits exhibited by insect societies 
seem closely related to phase transitions. A wide range of examples 
have been analyzed, from worker specialization (Holldobler and 
Wilson 1990) to nest architecture (Valverde et al. 2009). In the 
latter case, it has been shown that the spatial structures associated 
with the organization of termite nests seem to be compatible with 
a percolation transition. 


15 


LANGUAGE 


15.1 Lexical Change 


Human language stands as one of the most important leaps in 
evolution (Bickerton 1990; Deacon 1997; Maynard Smith and 
Szathmary 1995). It is one of evolution’s most recent inventions, 
and might have appeared as recently as fifty thousand years ago. 
Our society emerges, to a large extent, from the cultural evolution 
allowed by our symbolic minds. Words constitute the substrate 
of our communication systems, and the combinatorial nature 
of language’ (with an infinite universe of sentences) allows us 
to describe and eventually manipulate our world. By means of 
a fully developed communication system, human societies have 
been able to store astronomic amounts of information far beyond 
the limits imposed by purely biological constraints. As individuals 
sharing our knowledge and the cumulative experience of past 
generations, we are able to forecast the future and adapt in ways 
that only cultural evolution can permit. 

The faculty of language makes us different from any other 
species (Hauser et al. 2002). The differences between animal 


! More important, language involves recursivity, by which very complex, 
coherent sentences can be constructed in a nested manner. 
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Figure 15.1. (a) Bifurcations in word learning dynamics. If the rate 
of word learning exceeds one (i. e. R; > 1), a stable fraction of 
the population will use it. If not, then a well-defined threshold is 
found leading to word extinction. The inset shows an example of the 
logistic (S-shaped) growth curve for R; = 1.5 and x;(0) = 0.01. 
Lexical diffusion also occurs among artificial agents (b) where words 
are generated, communicated and eventually shared by artificial, 
embodied agents such as robots (picture courtesy of Luc Steels). 


communication and human language are fundamental, both in 
their structure and function. Although evolutionary precursors 
exist, it is remarkable to see that there seems to be no intermediate 
stage between them (Ujhelyi 1996). Such a shift might have 
occurred for a number of reasons (Wray 2002) from rare events 
(making human language a rather unique, unlikely phenomenon) 
to so-called macromutations. But alternative scenarios support 
the idea that language evolved from simple to complex commu- 
nication involving error limits (Nowak and Krakauer 1999) or 
phase transitions (Ferrer and Solé 2003). 

Far from being a static structure, language is always changing 
(McWorter 2003). Even at short time scales of only hundreds of 
years, as can be seen from the study of written texts (figure 15.1a), 
languages change fast. Languages themselves appear and disappear 
(Crystal 2000) and in some ways they behave as living species. 
As it happens with the fossil record of life, most languages of 


15.2 Words: Birth and Death 169 


the world went extinct a long time ago. As John McWorter 
points out: “Like animals and plants, languages change, split into 
subvarieties, hybridize, revivify, evolve functionless features and 
can be even genetically altered.” Perhaps not surprisingly, these 
similarities suggest that it makes sense to think of language and 
its changes in ecological and evolutionary terms. Many different 
aspects of language change have been addressed by theoreti- 
cians, and evidence indicates that phase transition phenomena 
might actually contribute to how languages are shaped through 
time. Here we consider such transitions in two very different 
contexts. 


15.2 Words: Birth and Death 


The potential set of words used by a community is listed in 
dictionaries (Miller 1991), which capture a given time snapshot 
of the available vocabulary; but in reality speakers use only part 
of the possible words: many are technical and thus only used 
by a given group, and many are seldom used. Many words are 
actually extinct, since no one is using them. On the other hand, 
it is also true that dictionaries do not include all words used by 
a community and also that new words are likely to be created 
constantly within populations, their origins sometimes recorded. 
Many such words demonstrate new uses of previous words or 
recombinations, and sometimes they come from technology. One 
of the challenges of current theories of language dynamics is 
understanding how words originate, change, and spread within 
and between populations, eventually becoming fixed or extinct. 
In this context, the appearance of a new word has been compared 
to a mutation. 

As occurs with mutational events in standard population 
genetics, new words or sounds can disappear, randomly fluctuate, 
become or fixed. In this context, the idea that words, grammatical 
constructions, and sounds can spread through a given population 
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was originally formulated by William Wang to propose and 
explain how lexical diffusion (i.e., the spread across the lexicon) 
occurs. Such a process requires the diffusion of the innovation 
from speaker to speaker (Wang and Minett 2005). 

A very first modeling approximation to lexical diffusion in 
populations should account for the spread of words as a con- 
sequence of learning processes (Shen 1997; Wang et al. 2004; 
Wang and Minett 2005). Such a model should establish the 
conditions favoring word fixation. As a first approximation, let us 
assume that each item is incorporated independently (Shen 1997; 
Nowak et al. 1999). If x; indicates the fraction of the population 
knowing the word W;, the population dynamics of such a word 
reads: 


dx; 
dt 


= Rx; = xj) — Xi (15.1) 


with 7 = 1,..., 7. The first term on the right-hand side of the 
this equation introduces the way words are learned. The second 
deals with deaths of individuals at a fixed rate (here normalized 
to one). The way words are learned involves a nonlinear term 
where the interactions between those individuals knowing W; (a 
fraction x;) and those ignoring it (a fraction 1 — x;) are present. 
The parameter R; introduces the rate at which learning takes 
place. 
Two possible equilibrium points are allowed, obtained from 
dx;/dt = 0. The first is x* = 0 and the second 
1 
+ 
k == R; (15.2) 
The first corresponds to the extinction of W; (or its inability 
to propagate) whereas the second involves a stable population 
knowing W;. The larger the value of R;, the higher the number 
of individuals using the word. We can see that for a word to be 
maintained in the population lexicon, the following inequality 
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This means that there is a threshold in the rate of word prop- 
agation to sustain a stable population. By displaying the stable 
population x* against R; (figure 15.1b), we observe a well- 
defined phase transition phenomenon: a sharp change occurs at 
R; = 1, the critical point separating the two possible phases. The 
subcritical phase R; < 1 will inevitably lead to loss of the word. 

An important implication of this “epidemic” of word spread- 
ing is that words used below threshold will disappear. Their 
extinction might be slow, starting from their disappearance from 
the social context and surviving in dictionaries, only to eventually 
be removed altogether.” 


15.3 String Models 


A step further in understanding language dynamics would make 
use of lists of words. This is obviously an oversimplification, since 
language is not just a collection of objects but also a complex set 
of rules relating and combining them (Solé et al. 2010b). Here 
we ignore the fundamental role of grammar and reduce a given 
language to a population of users sharing the same set of words. 

A fruitful toy model of language change is provided by the 
string approximation (Stauffer et al. 2006; Zanette 2008). In this 
approach, each language £; is treated as a binary string, that is, 
Li = (Si, Si, ..., SŻ) of length Z. Here Si €e {0,1}, and, as 
defined, a finite but very large set of potential languages exist. 
Specifically, a set of languages £ is defined, namely 


L= [Li L£L2,..., £4} (15.4) 


? For examples, see http://www.savethewords.org. 
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Figure 15.2. String language model. Here a given set of elements 
defines a language. Each (possible) language is defined by a string of 
v bits (here L = 3) and thus 2} possible languages are present in the 
hypercube. The two types of elements are indicated as filled (1) and 
gray (0) circles, respectively. 


with M = 2”. These languages can be located as the vertices of a 
hypercube, as shown in figure 15.2 for L = 3. Nodes (languages) 
are linked through arrows (in both directions) indicating that two 
connected languages differ in a single bit. This is a very small 
system. As L increases, a combinatorial explosion of potential 
strings takes place. 

A given language £; is shared by a population of speakers, to 
be indicated as x;, and such that the total population of speakers 
using any language is normalized (i.e, )>,x;=1). A mean 
field model for this class of description has been proposed by 
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Damian Zanette, using a number of simplifications that allow 
an understanding of the qualitative behavior of competing and 
mutating languages (Zanette 2008). A few basic assumptions are 
made in order to construct the model. First, a simple fitness 
function (x) is defined. This function measures the likelihood 
of abandoning a language. This is a decreasing function of x, 
such that @(0) = 1 and $(1) =0. Different choices are possible, 
including for example 1 — x, 1 — x? or (1 — x)?. On the other 
hand, mutations are also included: a given language can change if 
individuals modify some of their bits. 

The mean field model considers the time evolution of 
populations assuming no spatial interactions. If we indicate 
x= (x1,..., Xm) the basic equations will be described in terms 
of two components: 


dx; 
dt 


= A;(x) — M(x) (15.5) 


where both language abandonment A;(x) and mutation M; (x) are 
introduced. The first term is chosen as: 


A; (x) = px; ($) — @:)) (15.6) 


for the population dynamics of change due to abandonment. This 
is a replicator equation, where the speed of growth is defined by 
the difference between average fitness (p), namely 


N 
ORDD AZ (15.7) 
j=l 


and the actual fitness @(x;) of the i-language. Here p is the 
recruitment rate (assumed equal to all languages). What this fitness 
function introduces is a multiplicative effect: the more speakers 
use a given language, the more likely it is that they will continue 
to use it and that others will join the same group. Conversely, if 
a given language is rare, its speakers might easily shift to some 
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other, more common language. 
The second term includes all possible flows between “neigh- 
boring” languages. It is defined as: 


N N 
LL 
Mix) = > |) Wixi xi Wii (15.8) 
j=l j=l 


In this sum, we introduce the transition rates Wj; of mutating 
from language £; to language L; and vice versa. Only single 
mutations are allowed, and thus W;; = 1 if the Hamming distance 
D(£;, L;) is exactly one. More precisely, if 

L F 

DL: L) =X I$- Sf] =1 (15.9) 

k=l 
then only nearest-neighbor movements through the hypercube 
are allowed. In summary, A(x) provides a description of com- 
petitive interactions whereas M(x) gives the contribution of small 
changes in the string composition. The background “mutation” 
rate u is weighted by the matrix coefficients Wj; associated with 
the likelihood that each specific change will occur. 

This model is a general description of the bit string approxima- 
tion to language dynamics. However, the general solution cannot 
be found and we need to analyze simpler cases. An example is 
provided in the next section. Although the simplifications are 
strong, numerical models with more relaxed assumptions seem 
to confirm the basic results reported below. 


15.4 Symmetric Model 


A solvable limit case with obvious interest to our discussion 
considers a population where a single language has a population 
x whereas all others have a small, identical size, that is, x; = 
(1—x)/(N— 1). The main objective for defining such a symmetric 
model is to make the previous system of equations collapse into 
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a single differential equation, which we can then analyze. In 

particular, we want to determine when the x = 0 state will be 

observed, meaning that no single dominant language is stable. 
As defined, we have the normalization condition, now de- 


fined by: 


N M-1 — 
n =] 15.10 
Ds Ce) (15.10) 


(where we choose x to be the M-th population, without loss of 
generality). In this case the average fitness reads: 


sao a eee 
(@) = 00) + aa) (G4) : 


j=l 


Using the special linear case @(x) = 1 — x, we have: 


A(x) = px(1 — x) (« = yS) (15.12) 


The second term is easy to obtain: since x has (as any other 
language) exactly Z nearest neighbors, and given the symmetry 
of our system, we have: 


B(x) = = (2 a - s2) =-y (=) (15.13) 


And the final equation for x is thus, for the large- limit (i.e., 
when N > 1): 
dx 


— = px*(1— x) — ux (15.14) 
dt 


This equation describes an interesting scenario where growth 
is not logistic, as happened with our previous model of word 
propagation. As we can see, the first term on the right-hand 
side involves a quadratic component, indicating a self-reinforcing 
phenomenon. This type of model is typical of systems exhibiting 
cooperative interactions, and an important characteristic is its 
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hyperbolic dynamics: instead of an exponential-like approxima- 
tion to the equilibrium state, a very fast approach occurs. 

The model has three equilibrium points: (a) the extinction 
state, x* = 0 where the large language disappears; (b) two fixed 
points x} defined as: 


(15.15) 
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As we can see, these two fixed points exist provided that 
U < H: =p/4. Since three fixed points coexist in this domain 
of parameter space, and the trivial one (x* = 0) is stable, the 
other two points, namely x* and x7}, must be unstable and stable, 
respectively. If u < me, the upper branch x¥ corresponding to a 
monolingual solution, is stable. 

In figure 15.3a we illustrate these results by means of the 


bifurcation diagram using p = 1 and different values of u. In 
terms of the potential function we have: 
x? x4 x2 
Oana a 15.16 
u(x) p(5 +5) (15.16) 


In figure 15.3a—d three examples of this potential are shown, 
where we can see that the location of the equilibrium point is 
shifted from the monolanguage state to the diverse state as ju is 
tuned. The corresponding phases in the (p, u) parameter space 
are shown in figure 15.3b. 


15.5 Language Thresholds 


It is interesting to see that this model and its phase transition is 
somewhat connected to the error threshold problem associated to 
the dynamics of RNA viruses. In order for a single language to 
maintain its dominant position, it must be efficient in recruiting 
and holding speakers. But it also needs to keep heterogeneity (re- 
sulting from “mutations”) at a reasonably low level. If changes go 
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Figure 15.3. Phase transitions as bifurcations in Zanette’s mean 
field model of language competition. In (a) we show the bifurcation 
diagram using u/p as the control parameter. Once we cross the 
critical point, a sharp transition occurs from monolanguage to lan- 
guage diversity. This transition can be visualized using the potential 
function ®,,(x) whose minima correspond to possible equilibrium 
points. Here we use p = 1 with (b) u = 0.1, (c) u = 0.2, and (d) 
y= 0:3: 


beyond a given threshold, there is a runaway effect that eventually 
pushes the system into a variety of coexisting sublanguages. An 
error threshold is thus at work, but in this case the transition is 
of first order. This result would indicate that, provided a source 
of change is active and beyond the threshold, the emergence of 
multiple unintelligible tongues would be expected. 

String models of this type only capture one layer of word 
complexity. Perhaps future models will consider ways of intro- 
ducing further internal layers of organization described in terms of 
superstrings. Such superstring models should be able to introduce 
semantics, phonology, and other key features that are considered 
relevant. An example of such is provided by models of the 
emergence of linguistic categories (Puglisi et al. 2008). 

The previous models are just a toy representation of language 
complexity. Human language is the outcome of both social 
and evolutionary pressures together with cognitive constraints 
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(Christiansen and Chater 2008). Brains and languages (as well 
as languages and societies) coevolve (Deacon 1997). In this 
sense, the models presented above deal with language dynamics 
without the explicit introduction of cognition. Although our 
picture can be useful in dealing with the ecology of language 
(Solé et al. 2010b), it fails to include the adaptive potential of 
communicating systems themselves. Moreover, languages cannot 
be encapsulated in terms of word inventories. Words interact 
within networks (Solé et al. 2010a) related to the correlations 
established among them (phonologic, syntactic, or semantic). 
Phase changes play a role here, too, from the emergence of word- 
word combinations (Nowak and Krakauer 1999) to language 
acquisition in children (Corominas-Murtra et al. 2009). 


16 


SOCIAL COLLAPSE 


16.1 Rise and Fall 


The reasons for the disintegration of social order in ancient 
civilizations have been a recurrent concern for archaeologists, his- 
torians, and economists alike (Tainter 1988; Sabloff 1997; Fagan 
2004; Diamond 2005; Turchin 2006). Social collapse involves 
some kind of (historically) sudden, major loss of economic and 
social complexity. The fall of the Roman Empire, the collapse of 
the Mesopotamian centers, and the loss of the Mayan civilization 
are well-known examples of such large-scale decay (Tainter 1988; 
Yofee and Cowgill 1988). In all these examples, collapse implies 
a transition toward a simpler, less organized society that can 
even disintegrate and vanish. Information flows and economic 
transactions among different subsystems decrease or even stop, 
and social coherence shrinks. 

The laws that drive social collapse (if such general laws exist) 
are not known yet, but many clues can be recovered from the 
analysis of the past. An example is provided by the study of 
the Chaco Canyon site, in the northwest corner of New Mexico 
(Axtell et al. 2002; Kohler et al. 2005). This was the land of the 
Anasazi, and is today a semiarid ecosystem of undeniable beauty, 
but a hard place to live. The site contains the ruins of several sets 
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Figure 16.1. (a) An example of the large stone buildings of Chaco 
Canyon (picture by Bob Adams, see http://en.wikipedia.org/wiki/ 
PuebloBonito). It shows the basic organization of the Great House of 
Pueblo Bonito and involves hundreds of rooms (to be inhabited and 
for other purposes). In (b) the fluctuations in Anasazi populations 
in Long House Valley (Arizona). Population estimates are based on 
the number of house sites at a given time, assuming five people per 
house. Broad changes are observed, with a trend of increase marked 
by rapid shifts and a catastrophic decay between years 1250 and 
1300, when the Anasazi completely abandoned the valley (redrawn 
from Axtell et al. 2002). 


of buildings (figure 16.1a) connected to other distant locations 
through a network of roads spanning many miles. 

Until its disappearance around 1350, the Anasazi society was 
marked by fluctuations in total population size, mean settlement 
size, and preferred habitat. It is generally agreed that the reasons 
for such fluctuations (as it happens with the Mayan case) depend 
partly on rainfall and on groundwater levels, variations in either 
of which would affect Anasazi maize-based agriculture. Indirect 
measures of both variables have been obtained for some areas in 
northeastern Arizona (Figure 16.1b), together with archaeological 
surveys of the numbers and distribution of houses for the entire 
span of the valley’s Anasazi occupation. But how can one assess 
whether the proxies explain the survey results? To give just a 
hint of the complexities involved: readily available groundwater 
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can make rainfall unnecessary for maize agriculture, but the 
dependence of maize growth on rainfall and groundwater levels 
varies among the valley’s six habitats used for agriculture. 

It might sound like a bold idea trying to compress historical 
dynamics of any kind into a single equation of the form dx /dt = 
f(x). But properly formulated models and well-chosen variables 
have been shown to capture the key features of social and ecolog- 
ical dynamics (Levin 1999; Diamond 2002; Turchin 2003). The 
next example provides a good illustration of this point. 


16.2 Slow Responses and Sunk-Cost Effects 


The model presented here is based on the work described in 
(Janssen et al. 2003). This study provided a rationale for the 
presence of both resource depletion and failure to adapt in ancient 
societies. These models are based on the presence of a well- 
known deviation from rational decision making. Such a scenario 
is known as the sunk-cost effect (Arkes and Ayton 1999). The key 
ingredient in this model approach is an important feature of this 
effect: the fact that political decisions in social groups are achieved 
by some form of consensus. But once consensus is reached, the 
group will tend to commit to the decision taken, even if negative 
results are at work. 

Such a situation can generate an important and dangerous 
inertia. Quoting Jansen et al.: “even when the group is faced with 
negative results, members may not suggest abandoning an earlier 
course of action, since this might break the existing unanimity.” 

More generally, the underlying problem here is why complex 
societies might fail to adapt to resource depletion. Even if there is 
some social perception of risk, short-term thinking often prevails 
when facing long-term vulnerabilities. Such undesirable behavior 
is often favored by a combination of incomplete understanding of 
the problem, together with the misleading view that all changes 
are reversible. 
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Slow social responses to ongoing problems while in their early 
phases can result in a high overall cost (Scheffer et al. 2003). 
More precisely, lags in social response to a changing world can 
cause late, sudden shifts. Thus societies that respond slowly to 
challenges might be unable to avoid the long-term consequences. 
Delays between recognition and regulation can result from 
conflicts among individuals having different perceptions, social 
positions, and wealth. The inertia that results from these internal 
conflicts can easily lead to delayed decision-making dynamics. 

The consequence of the scenario described above is that 
there is no adaptive response to the changing resource and as a 
consequence the global dynamics can be described in terms of 
a resource-consumer system where resources are exploited with 
no particular forecast of future changes. As will be shown below, 
collapse can be a natural outcome of such a trend. The main threat 
here is a growing human population blindly exploiting limited 
resources. 


16.3 Ecological Model of Social Collapse 


Let us consider a simple description of the sunk-cost effect as 
defined by a one-dimensional model. For our analysis we take the 
approach used in (Janssen et al. 2002). The equation considers 
only resources x as the key variable and human population H is a 
fixed (but tunable) parameter: 

dx 


x 
Z = bx (1 z =) -T(x)H +8 (16.1) 


where u and K are the growth rate and maximum level of the 
resource, H is the current human population exploiting this 
resource, and ô indicates some resupply of x from neighboring 
areas. 

The functional response T(x) should introduce the way 
resources are depleted by humans, who act as consumers having 
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a given per capita rate of consumption C. Janssen et al. use the 
following response: 
Cx 
p+x 


r(x) = (16.2) 
where p is the resource level needed to reach 50 percent of the 
maximum consumption rate. The fixed points of this dynamical 
system can be obtained, as discussed in chapter 2 by looking at 
the intersection of two curves: 


1- Ž)+8 ple) =TQ)H (16.3) 


Those x* such that y;(x*) = y2(x*) will be the equilibrium 


yilx) = ux ( 


points. 

For simplicity, let us consider the case 6 = 0 (i.e., no external 
resupply). Such a simplification allows us to easily find the fixed 
points. The first corresponds to the (trivial) state x* = 0 
representing the complete depletion of resources. The stability of 
this fixed point can be determined from the sign of: 


_ _ 2x B pCH 


We obtain à (0) = u — CH/p, which will be negative if the 
human population density exceeds a threshold: 


ph 


H>H=— 16.5 

> C (16.5) 
The other two points are obtained from 
x CH 

1-—])- =0 16.6 

H( 3 p+x* a 


which gives two possible solutions, namely 


i=} |e-p ar ot si (o- )] (16.7) 


x 
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Figure 16.2. Phase transition in the sunk-cost model. In (a) we show 
several curves corresponding to the transient dynamics of x(¢) for 
increasing values of H (growing from top to bottom). At a given 
critical value H,, the available resources collapse. The corresponding 
bifurcation diagram is shown in (b). Here we have used K = 6, 
p = 1, u = 2, and C = 1, which gives a bifurcation for H, = 4.08. 
For H < H = 2 the unstable branch does not exist and small values 
of x lead to a rapid flow towards high x*. 


* 


The two fixed points x4 will exist provided that the sum inside 


the square root is positive (otherwise it would give imaginary 
values). This leads to a critical value of H that defines the domain 
of existence of nonzero resources: 

gat + 0)" (16.8) 

4KC 

The fixed point xj. is associated with a large amount of resources 
and thus a sustainable human population. In figure 16.2a we 
show several examples of the trajectories followed by our system 
for both H > H., and H < H.. But what about the intermediate 
one, x*? It corresponds to an unstable fixed point, separating two 
alternative states, as indicated in figure 16.2b with a dashed line. 
We can see that the only fixed point for H < py/C is actually 
xj, thus indicating that a small amount of resource would allow 
a rapid recovery. 
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Figure 16.3. Transitions from sustainable exploitation to collapse. 
In (a) the potential functions associated to the sunk-cost model are 
shown, for three different values of the bifurcation parameter H 
(other parameters as in the previous figure). In (b) the two phases 
exhibited by the sunk-cost model. The critical line H, is represented 
against the per capita consumption rate C. The gray area corresponds 
to the collapse of available resources. 


The shifting dynamics between sustainability to collapse can 
be described by means of the potential ®(x), which in our case 
reads: 

U, H 


(x) = =" + rail —CH[x+pln(o+-x)] (16.9) 


This potential function is displayed in figure 16.3a for different 
values of human population densities H. As this parameter grows, 
there is a continuous change in the depth of the valley associated 
to the nonzero state. Once the threshold is crossed, there is a 
transition from a finite resource level to its collapse (figure 16.3a). 
The critical boundary (16.8) allows defining the two key phases 
for our system (figure 16.3b). These are sustainable growth and 
social collapse, respectively. Either increasing population size at 
a fixed consumption rate or introducing higher consumption 
at a given population level (i.e., more efficiency in resource 
exploitation) will eventually push the system into the collapse 
phase. Although we have assumed that no external income of 
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resources is present (i.e., 6 = 0) the basic conclusions reached 
here remain the same. As the reader can check, introducing 
increasing levels of resource income allows the lower branch to 
involve a finite amount of stable resources. 


16.4 Historical Dynamics 


The previous model needs to be taken (as all models ana- 
lyzed here) as a first approximation to reality. It is more likely 
that the patterns of collapse and the conditions for sustainable 
development are themselves nonlinear functions of different social 
features. The model considered here lumps together a wide range 
of economic and sociologic traits. It also reduces the potential 
diversity of available resources to a single, average resource com- 
partment. Nevertheless, the study of the past and the possibly 
unavoidable presence of rapid responses to continuous changes 
suggest that the basic message is robust. 

What might be missing here? One important component is 
the presence of extreme events associated with the fluctuations 
in external factors. This was the case when a series of severe 
droughts affected the Yucatan Peninsula and led to multiyear 
periods of very low precipitation. The Maya collapse seems cor- 
related with several waves of drought (Gill 2000). Coupled with 
population growth near the limits of available resources, environ- 
mental degradation and warfare combined to doom Mayan cities 
(Diamond 2005). Not a single reason seems enough to fully 
account for this outcome. But in general, a basic pattern can be 
delineated involving the previous ingredients plus social inertia. 

Collapse is only part of the historical dynamics of social and 
economic structures (Turchin 2006; Krugman 1996; Bouchaud 
and Cont 1998; Lansing 2005). The rules governing social 
dynamics have been a rich source of theoretical and statistical 
analysis. Statistical physics and agent-based approaches have been 
particularly useful in providing insight (Axelrod 1984; Epstein 
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and Axtell 1996; Weidlich 2000; Bettencourt et al. 2007; Ball 
2008). Many aspects of social organization can be measured and 
scrutinized by means of appropriate theoretical models. Examples 
include the distribution of wealth, city formation, public opinion, 
economic fluctuations, spatial segregation or war (Choi and 
Bowles 2007). Is it possible to formulate a quantitative theory of 
history? Some authors suggest that such dynamical theory could 
be a reality (Homer-Dixon 2006; Turchin 2007). In particular, an 
evolutionary view of cooperation and conflict among individuals 
suggests that the former helps building societies whereas the 
second (once societies become empires) favors their dissolution. If 
true, our view of history as a dynamical, large-scale phenomenon, 
would benefit from incorporating the concept of phase transitions 
and their associated warning signals. 
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168-69; evolutionary origin of, 
167-68; limitations of simple 
models for, 177-78; recursivity in, 
167n; word diffusion in, 168, 
169-71 

language acquisition, 10, 178 
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lattices: Bethe lattices, 59-61, 60; 
percolation on, 53-57, 54, 55, 
59-61, 60; spins on (see Ising 
model) 

learning, of words, 168, 170-71 

Leptothorax, 162, 163 

lexical diffusion, 168, 169-71 

life origins, 70-71, 75-77. See also 
molecular replicator model 

linear differential equations, 26-28; 
limitations of, 30; system of, 
29-30 

linear dynamical systems, 26-28 

linear stability analysis (LSA), 33-36 

links (edges), 63. See also degree of a 
node 

logistic growth model, 31-32, 32; in 
cancer instability analysis, 130-31; 
stability analysis of, 36; of word 
learning, 168 


magnetization in Ising model, 12, 
14-15, 15; critical exponent f 
and, 17; fluctuations in, 15, 
15-16; mean field equation for, 
21; symmetry breaking and, 39 

Malthusian dynamics. 

See exponential growth 

marginal state, 25, 26, 33, 35 

Mayan civilization, 179, 180, 186 

mean field models: advantages of, 23; 
Ising model in approximation of, 
18-23, 19, 22; stability theory 
and, 36 

measles, 103, 104, 107 

medicine, sudden transitions in, 136 

Mediterranean ecosystems, vegetation 
in, 137 


membrane receptors, 11 
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memory effects: broken symmetries 
and, 52; hysteresis and, 48; lake 
turbidity and, 135 

Mesopotamian social collapse, 179 

metabolism, random network models 
of, 64 

metapopulation: defined, 143; 
dynamics of, 143-46, 144, 146 

Metropolis algorithm, 14 

microtubules, 92, 92-98, 95, 97 

molecular interactions, and phases of 
matter, 6 

molecular replicator model, 71-75, 
73; 75; 76 

Monte Carlo simulation, of Ising 
model, 13-17, 15 

morphogenesis: sol-gel transition and, 
58-59. See also differentiation, 
cellular 

multicellular systems, 120-21 

mutation rates: in cancer cells, 89, 
121, 130; critical, of quasispecies 
model, 83, 84, 86-88; error 
threshold of, 80; vs. genome size, 
79, 87-88, 88; of RNA viruses, 
78-80, 79, 88, 89; robust genomes 
and, 90 

mutations: cancer associated with, 
120, 128; in languages, 169, 173, 
174, 176-77; normal mechanisms 
for controlling, 121, 128; 
selectively neutral, 89-90 

mutator phenotype, 128 


natural selection. See selection 

networks: heterogeneous, 68-69, 
107, 155—56; infectious diseases 
and, 107-8; neutral, of genotypes, 
89-90; origin of life and, 77; 
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random homogeneous, 63-68, 65, 
66, 68; real, 68—69, 107, 155-56; 
robustness of, 69; scale-free, 69, 
108; of words, 178. See also 
computer networks 

neural networks, 64, 117 

neutral networks, 89-90 

nodes (vertices), 63. See also degree of 
a node 

nonequilibrium phase transitions, 
23-24, 57, 61-62 

nonlinear differential equations, 
30-33 


nonlinear dynamical systems, 30-33 


ocean sediments, dust in, 44—46, 45 

Onsager, Lars, 19 

opinion formation, 10 

optimization, based on collective 
behavior, 159, 165 

order, emerging in random graph, 68 

order-disorder transitions, in ant 
colonies, 162-65, 163, 165 

order parameter, 11; average 
magnetization as, 14—15; 
continuously changing, 44; for 
percolation on lattice, 56 

ordered phase, of Ising model, 12 

oscillations. See fluctuations 


pandemic, 108 
parallel multiprocessor networks, 149 
parasites, 78; of extinct species, 136 
partition function, of Ising model, 13 
path-dependent phenomena, 38 
patterns: complexity and, 2; 
symmetry breaking and, 38. 
See also spatial structures 
percolation, 53—62; as absorbing 


phase transition, 57; on Bethe 
lattices, 59-61, 60; bifurcation 
analysis compared to, 61-62; 
contexts of, 57; defined, 55; of 
neutral genotypes in sequence 
space, 89; in random Boolean 
networks, 113, 174, 115-17, 117; 
on random graphs, 63, 65-68, 68; 
as second-order phase transition, 
56; sol-gel transition and, 57-59, 
58, 98; on square lattices, 53-57, 
54, 55; termite nest organization 
and, 166 

percolation threshold, 63; for 
autocatalytic sets, 77; on Bethe 
lattices, 60; divergence close to, 61; 
for epidemic state, 100; in random 
Boolean networks, 116; on 
random graphs, 63, 65-68, 68; on 
square lattices, 55—56 

phase diagrams: mean field approach 
and, 23; nonequilibrium systems 
and, 24; for phases of matter, 4-7, 
5, 7; for water, 5, 6, 7 

phase transitions: absorbing, 57, 
100n; bifurcation analysis and, 38, 
61-62; examples in various 
contexts, 7—11, 8; in history, 187; 
in insect societies, 158, 159, 161, 
165, 166; language change and, 
168, 169, 171, 176, 177, 178; 
nonequilibrium, 23-24, 57, 
61-62; origin of life and, 77; in 
physics, 2-3; in tumor growth, 
120. See also critical boundary; 
critical point; first-order phase 
transitions; mean field models; 
networks; percolation; 
second-order phase transitions 
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phases: bifurcation diagram and, 41; 
defined, 2 

phases of matter, three standard, 2-3, 
3; phase diagrams for, 4-7, 5, 7 

pheromones, and foraging by ants, 
159-61, 160, 162 

Poincaré, Henri, 37 

Poisson distribution, 65, 67 

politics, 10 

polymers, 58, 58-59, 91, 92n, 
97-98. See also microtubules; 
proteins 

population growth: Allee effect and, 
145-46, 146; resource depletion 
and (see sunk-cost effect) 

population shifts, sudden, 134 

potential function, 41-43, 43 

power grid, congestion in, 148 

power laws, 16—17; for clusters in 
Ising model, 16; critical mutation 
rate and, 88, 88; real networks 
and, 69 

prebiotic evolution, 70-71, 75-77. 
See also molecular replicator 
model 

probability of a spin state, 13 

proteins, 58, 59, 91, 92n, 97-98; 
folding of, 8, 9; gene regulation 
and, 118-19, 119 


quasispecies, 80; bit string model of, 
81-83, 82, 83, 84; further 
developments of, 89-90 


radioactive decay, 28, 28n 

rainforests: ants and termites in, 157; 
fires in, 147 

Random Boolean Networks (RBN), 
112-13, 114, 115-17, 117 


random graphs, 63-69, 65, 66, 68 

recombination, and fitness 
landscapes, 89 

relaxation time, close to criticality, 
43, 49-51, 50, 165 

replicator equation, 173 

replicators, molecular, 70, 71-75, 73, 
75, 76 

resource depletion in societies, 
181-86, 184, 185 

RNA viruses, 78—80, 79; mutation 
rates of, 78-80, 79, 88, 89; neutral 
mutations in, 89—90 

robots, lexical diffusion among, 168 

Roman Empire, fall of, 179 

routing tables, depth of, 155-56, 
156n 


scale-free networks, 69; infectious 
diseases spreading in, 108 

scale invariance, 17 

second-order phase transitions, 
43-44; to epidemic state, 100, 
101; in Garay-Lefever cancer 
model, 123; in green-desert 
transition, 139, 139; percolation 
as, 56; of random network, 67. 
See also phase transitions 

selection: cancer cell population and, 
131; of prebiotic molecules, 71; in 
Swetina-Schuster quasispecies 
model, 86; viral mutation and, 88, 
89-90. See also evolution 

self-organization: of computer 
networks, 154-55; emergent 
phenomena and, 2; of insect 
societies, 157-58, 159; of 
vegetation in space, 136-37, 136n, 
137 
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self-organized criticality, 154-55 

semiarid ecosystems, 45, 134-35, 
137 

sexual contact networks, 108 

SIS model, 102-5, 104, 105 

smallpox, 99, 107 

social collapse, 179-86; factors 
predisposing to, 186; historical 
examples of, 179-81, 180, 186; 
resource depletion and, 181-86, 
184, 185 

social insect behavior, 157—66; 
broken symmetry in, 158, 158-61, 
160, 162; order-disorder 
transitions in colonies and, 
162-65, 163, 165; overview of 
complexity in, 157-58, 165-66 

social interaction networks, 107, 108 

social systems: catastrophes in, 44; 
Ising model of, 11; modeling of, 
187 

sol-gel transition, 57-59, 58, 98 

spatial correlations: in ecosystems, 
147; in epidemic spreading, 
101-2; ignoring, 23, 63, 68 
(see also mean field models) 

spatial dynamics: fitness landscapes 
and, 89; of infectious diseases, 107; 
of tumor growth, 120, 133 

spatial segregation, in human 
societies, 187 

spatial structures: broken symmetries 
and, 52; percolation and, 166; of 
vegetation, 136-37, 136n, 137. 
See also patterns 

spins, in Ising model, 12-13 

stability, 25, 26, 32-36 

stable equilibrium points, 25, 26, 33, 
34; of logistic model, 36; as 


minima of potential, 42, 43. 
See also fixed points 

stem-cell switches, 119 

sunk-cost effect, 181-82; model 
based on, 182-86, 184, 185 

superstring models, 177 

surfaces, growing, 10 

survival of the flattest, 90 

sustainable development, 184, 
184-86, 185 

swarm behavior: in fish schools, 8, 
9-10; in insects (see social insect 
behavior). See also cooperation 

Swetina-Schuster quasispecies model, 
84-88, 85, 88 

symmetry, and phases of matter, 
2-3 

symmetry breaking, 38-43, 39, 41; in 
ant colonies, 158, 158-61, 160, 
162; applications of, 52; critical 
slowing down and, 49-51, 50; 
fluctuations close to, 43, 44; in 
Ising model, 39, 41; multiple steps 
of, 51-52, 52. See also bifurcations 

system: defined, 53; threshold for 
existence of, 56 (see also percolation 
threshold) 

system of linear differential equations, 
29-30 

systems biology, 111 


Taylor expansion, 34, 34n 

technological evolution, 38 

temperature: as control parameter, 3, 
11, 15; critical, of Ising model, 15, 
15, 16, 17-18, 19; phases of 
matter and, 2—3 

tensegrity, 97—98 

termites, 158, 166 
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thresholds: bifurcations and, 38. 
See also eradication threshold; error 
threshold; percolation threshold 

traffic. See computer networks 

transient time, 49-51, 50; for 
percolation on lattice, 57 

transportation networks, 107-8 

trees, on square lattice, 53-57, 
54, 55 

tree structure, of Bethe lattice, 
59-60, 60 

tristable dynamics, 119 

tubulin, 92, 93, 94 

tuning of a parameter: bifurcation 
and, 38, 41. See also control 
parameter 


Turing-like instability, 136n 


universality, 17—18 

unstable equilibrium points, 25, 26, 
33, 34-35; of logistic model, 36; as 
maxima of potential, 42, 43. 
See also fixed points 
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vaccination, 105-7, 106 

vegetation patterns, self-organized, 
136-37, 136n, 137. See also 
green-desert transitions 

vertices (nodes), 63. See also degree of 
a node 

viruses, 78-80, 79; digital model of, 
81-83, 82, 83, 84; HIV, 79, 89; 
spreading in networks, 108. See 
also epidemic spreading; 
quasispecies; RNA viruses 

von Neumann neighborhood, 54 


Waddington, Conrad Hal, 51, 52 

war, 187 

water: amphiphilic molecules in, 9; 
first-order phase transitions in, 11, 
44; phase diagram for, 5, 6, 7 

water in ecosystems: Anasazi society 
and, 180-81; availability of, 134, 
135, 136 (see also green-desert 
transitions); Maya collapse and, 


186; turbidity shift in lakes, 135 


